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ABSTRACT 


When the thermal regime of permafrost is altered,melting of 
the frozen ground may be induced. Many civil engineering problems are 
concerned with the stability and deformation of such thawing soils. 
Some available solutions to the heat transfer problem are reviewed, 
which facilitate the calculation of the rate of melting. The physics 
of consolidation for thawing soils are formulated in terms of known 
equations for heat transfer and for linear consolidation. A moving 
boundary consolidation problem results, and analytical solutions are 
obtained for several cases of practical interest. The results are 
presented as normalised pore pressure distributions and settlement 
ratios, and it is demonstrated that both the pore pressures and 
settlements depend primarily on the thaw-consolidation ratio R. 

An apparatus was fabricated to simulate the required stress 
and thermal conditions. in the thawing soil. Samples of natural and 
reconstituted frozen soils were tested to determine the validity of 
the theoretical predictions, and extremely encouraging results were 
obtained. The concept of the residual stress is introduced, and its 
importance in thaw-consolidation and the undrained thaw strength 
is discussed. 

The versatility of the theory is greatly increased by 
considering several theoretical extensions which permit the analysis 


of a wide range of practical problems. In particular, solutions for 
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arbitrary movement of the thaw plane are obtained numerically, and 
some non-linear constitutive relationships for the soil skeleton are 
incorporated in the theory. The use of vertical sand drains for 
improving the conditions in a thawing foundation is discussed, and a 
typical case is analysed. 

The excess pore pressures, settlements and rate of thaw in 
the foundation of a hot oil pipeline at Inuvik, N.W.T., are analysed 
within the context of the thaw-consolidation theory. Excellent 
agreement is obtained between theoretical predictions and the 


observed behaviour in this case study. 
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CHAPTER I 


INTRODUCTION 


1.1 Permafrost and Associated Geotechnical Problems. 


Approximately one half of Canada is underlain by permafrost. 
Permafrost is more correctly termed "perennially frozen ground", and 
by definition refers to ground that remains frozen for more than one 
complete year. The thickness of permafrost may vary from a few feet 
at the southern boundary to greater than two thousand feet in the 
far north. The extent and many of the surface features associated 
with permafrost are described in detail by Brown (1970). 

The delicate thermal equilibrium of permafrost is particularly 
sensitive to any natural or man-made changes in envircnmental conditions. 
The rate of construction of buildings, transportation and storage 
facilities necessary for the development of the north and the 
exploitation of its resources is increasing continually. 

The problems in soil mechanics likely to be encountered in 
the Arctic are usefully subdivided into two sections, those involving 
a change of state, and those where no change of state takes place. The 
latter class of problems involves soil which remains continuously 
frozen or thawed. For geotechnical problems involving thawed soil, 
the conventional and well-established principles of soil mechanics 
are applicable. 


If the region of soil under consideration is continuously 
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frozen, the principal engineering problems are involved with the significant 
time-dependent deformation of frozen soils under constant stress. An 
engineering theory of creep for frozen soils has been established by 

Vialov (1963) and by Ladanyi (1972), based on a collection of 

rheological laws which are found to describe adequately the observed 
manifestations of creep at least in the laboratory. The theory as 

it presently stands seems perfectly capable of solving some specific 

soil engineering problems such as the bearing capacity of buried 

footings and anchors. 

More dramatic problems arise when the geotechnical engineer 
is confronted with a soil whose water phase is actively freezing or 
thawing. Examples of soil freezing may be seen in frost penetration 
under highways, the long-term aggradation of permafrost, the annual 
freeze-back of the active layer, and in the freezing of soil around 
chilled pipes or cold gas pipelines. It is known that most saturated 
soils under low stresses exhibit an affinity for water during 
freezing, and the formation of lenses and layers of ice may occur. 
Coarse-grained soils under large stresses may also expel water on 
freezing if drainage is permitted (Mackay, 1971). If ice lensing jis 
inhibited, substantial heaving pressures can result, or conversely, 
if the overburden stress at the frost plane is high enough, the 
growth of segregated ice is impeded. In reviewing the geotechnical 
properties of freezing ground, Anderson and Morgenstern (1973) 
have stated that ice lensing results when a soil has the ability to 
supply water to an active ice front for a sufficiently long time to 


grow an ice lens. The particle size of the soil influences the water 
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flow to the ice front partly through its effect on the permeability 

of the unfrozen soil, and perhaps more importantly, by influencing 

the maximum stress difference that can be Supported at the ice-water 
interface. This maximum stress difference is defined as the difference 
between the total and pore water stresses at the freezing plane, and 
has been the subject of much research, e.g. Williams (1966), Penner 
(1967). 

The problems involving freezing soils might present themselves 
in two ways. If ice lensing occurs at low stress levels, and remains 
unrestrained, then serious deformations may occur in structures such 
as highweys or chilled gas pipelines. On thawing of fine-grained soils 
containing segregated ice structures, serious stability and settlement 
problems may arise. If on the other hand, the ice lens growth is 
restrained during freezing, severe heaving pressures may arise which 
could cause damage to engineering structures. 

Finally, of the important problems confronting the 
geotechnical engineer in the Arctic at the present time, probably the 
most common and potentially serious are those involving the thawing 
of frozen soil. Adverse conditions may disturb the delicate thermal 
equilibrium of permafrost sufficiently to induce melting. Examples 
of this are the burning or stripping of vegetation, the impounding 
of a reservoir, the construction of heated buildings or the thermal 
erosion adjacent to northern rivers. An extreme case of the melting 
of permafrost has been shown to exist around a hot oil pipeline 
buried in permafrost (Lachenbruch, 1970). 


Severe conditions might be anticipated in a thawing 
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foundation if the soil is rich in ice, or the soil is fine-grained and 


the rate of degradation is rapid. 


1.2 Scope of the Thesis 

In more northerly locations, the engineer is principally con- 
cerned with ground that is initially frozen (permafrost) and which is 
in danger of thawing. Thus, an analysis involving the stability and 
deformation of thawing soil appears to have the highest priority for 
the geotechnical engineer who is working in an Arctic area. 

The analysis of deformation under stress, and the stability of 
the thawing Soil against shear failure are both known to be intimately 
related to the effective stress in the thawing soil. Once the effective 
stresses are known, combination with a constitutive relation for the 
soil yields the deformations, and combination with the conventional 
strength parameters provides the available Shear strength. 

To the present time, no correct or coherent theory is 
available to accomplish the analysis mentioned above for thawing soils. 
With the recent accelerating development in Arctic regions, a rational 
analysis for the behaviour of thawing Arctic soils takes on a high 
priority indeed. It is to this end that the research work described 


in this thesis was carried out. 


1.3 Review of the Geotechnical Literature Concerning Thawing Soils 


At the time of the First International Permafrost Conference 
(1963) almost no North American geotechnical studies had been published 
on the behaviour of permafrost subjected to thawing. The general 


characteristics of high compressibility and low shear strength had 
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been recognised in soils which were ice-rich and fine-grained prior 
to thawing. It was intuitively suspected that if the rate of generation 
of water in a thawing soil exceeded its drainage capacity, excess 
pore pressures would be generated, causing the reduction of the available 
snear resistance. This was the theme of a simplified analysis undertaken 
by Yao and Broms (1965), who recognised that the bearing capacity of 
a subgrade soil might be reduced by the incomplete dissipation of 
excess pore pressures during thawing. Their theoretical study did not, 
however, embrace fully the properties of fine-grained soils. An 
analysis for thawing and the subsequent total settlements under dikes 
at Kelsey, Manitoba was carried out by Brown and Johnston (1970). 
The total settlements were calculated on the basis of a visual 
estimation of the segregated ice conditions, but no analysis was 
undertaken for the transient settlements and excess pore pressure 
conditions. However, aS no serious stability problems arose subsequent 
to thawing under the dikes, it was concluded in retrospect that 
adequate drainage was available to Fate a stable foundation. 
In the North American literature, physical descriptions of 
the processes involved in the consolidation of thawing soils were 
presented by Aldrich and Paynter (1953) and by Scott (1969), but 
no rigorous statement of the problem or solutions were obtained. 
Considerable research has been continuing in the USSR 
Since about 1930. Although many papers of a descriptive nature were 
published (Tsytovich, 1957; Kiselev et al, 1965; Tsytovich, 1965) 
little analysis of the problems of thawing soils was ever attempted. 


The first analytical solution to a problem in thaw-consolidation 
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was apparently obtained by Tsytovich et al (1965). A situation was 
envisaged where the thaw front was penetrating into the soil at a 
rate proportional to the square root of time. The one-dimensional 
Terzaghi equation of consolidation was assumed to govern the 
dissipation of excess pore pressures. However, the analysis failed 
to describe the thaw-consolidation process in terms of previously 
defined geotechnical properties, because the excess pore pressure 
at the thaw line was assumed equal to a constant value. This constant 
pore pressure boundary condition was assumed to be equal to a fraction 
of the overburden stress on empirical grounds. This assumption 
necessarily required a laboratory test on thawing soil with pore 
pressure measurements, and was never related to other more easily 
recognised geotechnical parameters. This analysis, therefore, would 
be very difficult to use, and would probably only be applicable in 
practice to one set of soil, loading and thawing conditions. 

Another analysis for thawing soils was carried out by Feldman 
(1965), who attempted to incorporate a non-linear constitutive relation 
for the soil skeleton. Although admirable in intention, the analysis 
failed completely by assuming the effective stress at the thaw line 
was continually equal to zero. Furthermore, it was incorrectly 
assumed that all settlement took place during thawing. Consequently, 
this analysis adds little to the general understanding of the subject. 

A study by Malyshev (1966) incorporated a variation of total 
stress with depth, and although a sophisticated solution technique 
was required, no new physical statements were introduced to improve 


the behavioural description of thawing soils. 
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The first analysis which attempted to model realistically 
the boundary condition at the thaw line was proposed by Zaretskii 
(1968). Although the continuity of the pore phase at the thaw 
boundary was described mathematically in terms of recognisable soil 
parameters, errors were introduced in the formulation at this point. 
Even though the physical statement of the thaw front boundary 
condition was incorrect, this work provides the first valuable insight 
into the pore pressure and effective stress conditions in a thawing 
soil. 

Considerable interest in North America has recently been 
concerned with the proposed construction of a hot oil pipeline through 
Arctic regions. Such a proposal has intensified the research effort 
on geotechnical problems in Arctic Canada and Alaska. The thermal 
disturbance caused by a warm oil pipeline on the underlying 
permafrost was estimated by Lachenbruch (1970) in one of the most 
Significant studies carried out to date in this research area. This 
research established that a thawed area of some forty feet around the 
pipe was to be expected during its operating lifetime. Some of the 
implications were explored assuming that the thawed soil behaved 
like a viscous fluid, and it was demonstrated how this would be a 
serious impediment to the burial of a hot pipeline. This also 
provided a substantial incentive to develop a more realistic model 
for the behaviour of permafrost subjected to thawing. 

More recently, various agencies have constructed test pipeline 
sections, to determine the effects of a hot pipeline on the foundation 


material. The observed field behaviour of a test section at Inuvik, 


re 


eMtataatt en Fsbom ot porqmasts ote 2ieylens’ sits 
tetas al bszodn7g 26w ent “wens any ts ‘aotstbe Ey. 
werd eat os seéng stog aiif to’ itl sa0s ant 
fioe afdsetnpousy Yo 2ated ni yf ; sot oma an peer aa 
»tatog 2fdd te hore Tumvo?. add nk beoubontar SEM mond + 
“reboued fworxt wort ait to Jnamedste rte aha 
saptent seheniadd dantt art abt voNg Avow 2FHT ENT o8W - 
oatwallt 5 nf enottibnos ssovte evigosttg bing wiusoeud: 0g, 9 


: ve 
" ee 


need yidreoot 2ad sotvems rtdoll nf Jasvatat ssi 
Apuoids ‘eni-feqtq fro tod s 18 nat tourtetios bazoqorg ong a 
‘SOF, Horss297 943 borttenstar zed Feedaerg we dae 191 
fimuedt-ad isdasfA bas sbens) aftowé nb ei fdovy teat e 3 7 
ent Tvabnu itt no ontfearq: (Fo tirew 6 yal sanded cx | 

d2on Sat to ano nt (OXeT) aundnan ied yd: possniyes caw Jae ‘a 
2nat BONS donsseley 2tit ni etsb of +is0 bei ans 2atbude Faey he bi a 
aa bauows Hast (rer Smog to HIV bowed ® ead baieifdavas sy 
‘ots Yo Srm2 cemigarit pifiéydsao" ett boli betoagxe sd. ot | 
bauailed: Tree bawailt gut seis ghimees berolns avee eholtest he 
& a Bivow frit wor bosertencineh paw tt bas. belie es a 
dels atiT smi Feqig JON 6 6 fettid ont os {wonitbsn?’ ave ibe 
Faber otter tse vig 6 qofevab oti. sePnaont’ (stinatedue: é el 
_spatwedt od boy astade seoytoiyge 2 worveded ode all 

| an fagig s2ad BeSsuiseno. oved 2stsiSpe 2UOrnsy idneos4 oom 

motgebaue? att no Sni'fagty ton a 0 shoots efy aniaystss on ‘ela 


‘UR Wut 98 ebisobe esd '5 46 ‘woilae biger teva | a 


=) 


¢. 
“i 

> 

of 


N.W.T. was reported by Watson et al (1973). Settlements, pore 
pressures and thaw depths were reported, and these results constitute 
the first well documented field case history in North America on 
thawing permafrost that is available at this time. These observations 
are analysed in some detail in Chapter 6 of this thesis. 

However, research on some of the geotechnical properties 
of thawing soils is still limited on this continent. Speer et-al 
(1973) have reported Eieereculte of numerous thaw settlement tests 
on undisturbed permafrost samples exhibiting a wide variation of 
ground ice conditions. A statistical correlation was established 
between the total thaw settlement and the frozen bulk density of 
the material. The bulk density of the frozen soil appears to form 
a useful index to the quantity of total settlement in this instance, 
as it reflects the amount of ice and to a lesser extent the amount of 
air that is present in the soil pore spaces. A correlation such as 
this is useful in establishing a preliminary estimate of total 
settlement, but, of course, gives no indication of the transient 
pore pressure and settlement responses during thawing. 

As can be seen, there is an extreme paucity of information to 
date on the geotechnical behaviour of thawing soils. Moreover, there 
is not available at the present time a complete, correct and 
succinct statement of the physics of thawing soils. In addition, 
there is no available solution which might be applied to solve the 
engineering problems involving the deformation and strength of 
thawing fine-grained soils. This is surprising in the light of 


recent developments in Arctic areas, and provides considerable 
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stimulus to remedy the situation, at least in part. 
This thesis develops solutions to some of the more obvious 
problems in permafrost soils confronting geotechnical engineers 


at the present, and indicates the manner in which problems of a more 


complex nature might be solved. 


CHAPTER II 
THE ONE-DIMENSIONAL THAWING OF FROZEN SOILS 


2.1 The Thawing of a Homogeneous Frozen Soil Subjected to a Step 


Increase in Surface Temperature. 


An analysis of the engineering behaviour of a thawing soil 
requires a knowledge of not only the extent of the thawed soil, but also 
the rate of melting. The functional relationship between the depth of 
thaw and time is required before any geotechnical study of thawing soil 
is attempted. Once this information is extracted from the solution to 
the problem in heat transfer, additional information such as the 
temperature profiles are considered to be of second order importance. 

Some of the assumptions involving the uncoupling of the thermal and 
consolidation problems in this manner, and their validity are discussed 

in subsequent sections. In the following paragraphs, some classical 

solutions to the melting problem are reviewed for soils where a sudden constant 
increase in surface temperature has been applied. 

It is assumed that conduction is the chief mode of heat transfer 
within the soil mass. The mixture of soil particles and water is assumed 
to behave as a continuum, which is deemed valid considering that the 
scale of a soil particle is many orders of magnitude smaller than a 
representative mass of soil. Temperature boundary conditions are used 


exclusively in this study. Although it is quite feasible to include 
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boundary conditions of the heat flux type in a numerical study, this 
approach is still the subject of much active research, and beyond the 
scope of this thesis. 

As illustrated in Fig. 2.1, a uniform homogeneous frozen soi] 
is subjected to a step increase in temperature from Ty initially in the 
ground to Le at the surface. The properties of the frozen and thawed 
zone are considered to be homogeneous and independent of temperature. It 
is assumed further that the latent heat associated with the transformation 
from ice to water is liberated at a fixed melting point Te = 0. The 
analytical solution for the movement of the thaw front and the associated 
temperature fields has been found by Neumann about 1860. A statement of 
the equations and boundary conditions is given in concise form by Carslaw 
and Jaegar (1947), where the solution is also derived. The movement of 


the interface between the thawed and frozen soil zones is expressed by 
X=avt (2.1) 


where X is the depth of thaw 
t is the time 
and a is a constant which is determined as a root of the 


transcendental equation 
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INTERFACE TEMPERATURE = 0 


FROZEN SOIL 


Ke Cra Keene 


GROUND TEMPERATURE = Tg 


Fig. 2.1 The Neumann problem 
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K and Ke are the diffusivities of unfrozen and 
frozen soil [om?/s]. 

k and Ke are the thermal conductivities of unfrozen 
and frozen soil [ca1/°C cm sj], 

c. and Cp are the volumetric heat capacities of unfrozen 
and frozen soil [cal/°C cm? ], 

L is the volumetric latent heat of the soil [cal/em> soil}; 

T_ is the (uniform) initial ground temperature difference 
ce below freezing], 

iv is the applied constant surface temperature difference 

iC above freezing]. 
As the diffusivity of a zone is simply the conductivity divided 


by the volumetric heat capacity for that zone, then the solution for a 


appears to be a function of the seven variables: 


p> Cys Ces af Ts and L) (2.3) 


Qa = F(k k 
Their relative importance will be examined later. 


Equation (2.2) may be rewritten in such a way that only three 


dimensionless parameters emerge. 
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and Ste is the so-called Stefan number, and is defined as 

the ratio of sensible heat to latent heat by 

ey, 

L (2.6) 


Ste = 


The seven variables in eq. (2.3) are absorbed into three dimensionless 
numbers. This makes a graphical presentation of the solution to eq. (2.4) 
feasible. 

Equation (2.4) has been solved by the Newton-Raphson iteration 
scheme for finding the zeros of a function. It is of interest to note 
that the results are almost totally independent of the ratio VK) /Ke: 

The results for the normalized thaw rate, a/2v kK » are presented in 


Fig. 2.2 as a function of the two significant variables 
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Varying VK /Ke over a reasonably wide range of values does not produce 
a measurable difference in the normalized thaw rate. A reasonable value 
of VK [Ke for most soil conditions is 0.7, and this value was adopted 
for Fig. 2.2. This chart therefore represents a complete and accurate 
solution to the Neumann problem, with no assumptions other than those 
inherent in the original formulation. 

Considerable simplifications to the rigorous Neumann solution 


have been made by various authors. It is possible to set up a table 


oe | Bt 1 cee 
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of solutions in order of decreasing complexity and accuracy as in 

Table 2.1. This table applies only for the step increase in temperature 
at the surface. It appears that a formulation by Berggren (1943) has 
duplicated the earlier work of Neumann, and therefore is not discussed 
here. 

After the solution proposed with VK /Ke = 0.7, the next 
Simplification was made by Aldrich and Paynter (1953) in the so-called 
modified Berggren method, which assumes a slightly less realistic value 
of VK [Ke = 1.0. Moreover, slight errors are introduced by taking 


Come gC > which by definition implies also that Ke = Ke 
The solution may be greatly reduced in complexity if we assume 
that the temperature distribution in the frozen zone does not affect the 
rate of thaw. As will be seen later, this is a reasonable simplification 
ig: Ty is close to the melting temperature. Setting 1g = 0 in eq. (2.2), 


we obtain 


2 
Oo 
a 4k a 
aos e u erf = Ste (247) 
POV %e & =) 


Now on a semi-empirical basis eq. (2.7) may be approximated to 


a high order of accuracy by 


or} a Ste Ste 
py: ene? ie ee ( f sie) 
2/7 Ky (258) 


This equation is extremely simple to evaluate, and the results 
of eq. (Z2s/) and eq. (2.8) are plotted in Fig. 2.3. 

If a linear temperature distribution is assumed in the thawed 
zone, and the temperature profile in the frozen zone is again ignored, a 


solution may be obtained which is often used and is originally due 
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Or writing in dimensionless form 
Ont saty/ Ste. 
ave, ¢ 


The results from eq. (2.10) are also plotted for comparison 


(2.10) 


ine h LG eee 

The solution to the Neumann equation is presented in chart form 
in Fig. 2.2 and the o value can be extracted from it with a high degree 
of accuracy. If the ground temperatures are close to the melting point, 
then the simple function given by eq. (2.8) may be used to predict a, 
the relationship between depth of thaw and the square root of time. 

Having established these analytical expressions for the thawing 
of frozen soils, it is desirable to examine the thermal parameters 


involved and their calculation. 


2.2 The Thermal Properties of Soils. 


The simplifications introduced by various authors in Table 2.1] 
Suggest that some thermal properties may influence the thaw rate to a 
greater extent than others. For instance, the Stefan solution in 
eq. (2.9) depends only on kK) L, and (pe It is of interest therefore 
to assess the relative importance of the remaining four variables on 
the rate of thaw a. The Neumann solution will be used in combination 
with data for thermal conductivity for a saturated silt/clay soil taken 


from Kersten (1949). Now Ce, Cy, and L are also dependent on water 


(8.3) 
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20 
content, and they are derived from the simple relationships given later. 


So k ke, Co Ces and L can all be written as functions of the water 


Wes 
content w, and therefore the importance of each thermal property is 
assessed by plotting a against w, for a given step temperature T,. In 
examining each of these variables in turn, it has been assumed that 


the degree of saturation is 100%, and the specific gravity of soil 


SOltdsS 1S 2.7. 


(a) Ground Temperature, Ra 


The effect of the initial (uniform) ground temperature is 
determined over a range of values from 0 to -5°C. The effect of a 
change in ground temperature of habits £0 change the rate of thaw 
parameter o by about 2.5% (see Fig. 2.4). These calculations are based 
on an example where Te = +10°C. The effect of ground temperature on the 
rate of thaw can also be investigated by consulting the graphical 
Solution of the Neumann problem in Fig. 2.2. 

In problems where the ground temperatures are close to the 
melting point, it may often be reasonable to ignore them entirely. This 
will result in marginally overestimating a. 

(b) Thermal Conductivity of Frozen Soil, ke 

It is well known that the thermal conductivity of frozen soil, 
Kes is temperature dependent (Penner, 1970) and that this phenomenon 
is due to the presence of unfrozen water in fine grained soils below O°C: 


Because k._. is greater than k and, as there is only a slight 


ice 
temperature dependence of conductivity in soil minerals, Kersten (1949), 


water 


it can be seen that Kk. should be greater than ky FO SOLIS. Enact.) Tor: 
high water contents the ratio Ke/k,, should approach the ratio of 


conductivities of ice to water of 3.7. Kersten has also presented 
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relationships for ky and ke for silt/clay soils and his results are 
presented in Fig. 2.5. It can be seen that the ratio kK e/k, ranges 
mom bol tO. 37 

The effect of the absolute magnitude of ke on a may be studied 
by varying kK e/k from 1.0 to 3.7. At a given water content, the value of 
ky is used in the Neumann solution together with the ke value as defined 


by the ratio k ¢/k This Ke is taken to be constant over the range 0°Cc 


u 
to Ty The rather insignificant effect of absolute variations in Kk. 
on a is readily apparent in Fig. 2.6. 

It is suggested that for any thaw calculation the use of Ke = Kae 
= 5.3 mcal/°C cm s will yield adequately accurate solutions to a thaw 
problem except for the most highly plastic clays. It should be noted 
that Kersten's data on Ke given in Fig. 2.5 approaches Kine at a water 
content of 24% and then increases. This subsequent increase in Ke is 
felt to be incorrect and likely results from the function fitting method 
used. 

Problems may be studied by numerical techniques where ke is 
dependent on temperature and an example is presented in a later section. 
However, the effect of the temperature dependent ke on a can be studied 
using the analytical results already obtained. As has been shown in 
Fig. 2.6, varying ke from 1.0 ky LOO, Ky has little effect ona. The 
effects of any intermediate relationship kp = k, (6) on a will be bounded 
by the results already obtained for Ke = 1.0 ky and ke aes al Ka 

It may be concluded that neither the absolute magnitude nor the 
temperature dependence of the frozen conductivity ke has any significant 


influence on the determination of the rate of thaw a. The calculations 


involved in Fig. 2.4 and Fig. 2.6 were provided by McRoberts (1972). 
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(c) Volumetric Heat Capacity of Frozen Soil, Coda 


The volumetric heat capacity of soil is calculated by a simple 
mixture rule in which the relative heat capacity of the constituent 


materials are added together. That is 


cor * a (c+ 0-5 (2.11) 


where Yq denotes the dry density [g/cm?], 

CH is the heat capacity of the soil grains [cal/g.°C}, 

0.5 is the heat capacity of ice [cal/g.°C], 

w denotes the water content [g/g]. 

As pointed out by Moulton (1969), there are some minor differences 
of opinion among some authors as to the correct value of CH to use in 
eq. (2.11). Kersten (1949) finds that a value of 0.17 ca¥g.°C represents 
an average value for soils close to their freezing point. Moulton refers 
to several sources that suggest a value of 0.20 cavg.°C may be more 
representative. 

For frozen soils containing unfrozen water, a modification to 
eq. (2.11) can be made to allow for the relative percentage of ice and 


water as follows: 
Ce = Yq CH FuOLS (1 + W do (212) 


where Wh is the unfrozen water content [g. water/g. (ice + water)]. 
A value of Ce may be calculated based on the value of Wh at 

the ground temperature Lig of interest. However, as pointed out in the 

preceding sections, the effects of 1g and Ke on the rate of thaw a are 


small as well. For large ratios of T/T, the role of Ce in the overal] 
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solution can become more significant, but these ratios are generally 
not encountered in geotechnical problems. 
The effect of Ce can be inspected by considering the 


dimensionless variable 


es) ota a 
Li a Ke 


M1 1G. 92.2. AS Ke may be replaced by Ke/Ces it follows that a is related 
to ¥ Ce . For an example with w = 0.5 and We = 0.2 at lhe Ce is under- 
estimated by 10% if eq. (2.11) is used. This error results in an 
underestimation of the rane of thaw of less than 2% for the usual range 
Or “ote” ; 

If this error is considered to be significant, eq. (2.12) can be 
used to dering Ce. The effect of the unfrozen moisture content-temperature 
relationship between 1g and 0°C on Ce is completely insignificant with 


regard to predictions of the rate of thaw. 


(d) Volumetric Heat Capacity of Thawed Soil, Cree 


The parameter Cy plays an important role in thaw solutions. It is 
used in defining the Stefan number, Ste, which has been shown to be an 
important variable. The value of cj, can be determined, as was Ces by 


the following: 


Cay (, + 1.05) (Zk 3) 


where 1.0 = heat capacity of water [cal/g.°C}. 
The comments in the preceding section on Ce apply equally to Cu. 
Variations in CH have less effect on Cy than on Ces as the specific heat 


of water has doubled from that of ice. 
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(e) Thermal Conductivity of Thawed Soil, kK 


The rate of thaw a is proportional to the square root of ky 
and hence errors in ky have serious consequences on a. For example an 
werror, OT tb 25% in kK, results in an error in a of about + 11%. The 
relationship given by Kersten (1949) for the thermal conductivity of 


thawed soil is 


: =4 0.624, 
Ky 3.446 x 10 [0.9 10g, 9() os OeZ] 10 d (2.14) 


where ky is in cal/om.°C.s, 
w is the water content (%), 


Yq is the dry density (g/cm?) . 


Kersten's data on oF have been used for the parametric studies 
presented. It is felt, furthermore, that the absolute magnitudes of the 
ky values presented are entirely reasonable. Figure 2.5 indicates 
the comparison of Kersten's data with that presented by Makowski and 
Mochlinski (1956) and single determinations by Watzinger and Saare 


referenced by Skaven-Haug (1972). 


(f) Latent heat and Water Content Effects 

The water content of a soil is a dominant parameter in thermal 
calculations. As has been shown, the thermal conductivity and volumetric 
specific heats of both thawed and frozen saturated soil can be expressed 
as functions of the water content. Below 0°C it is recognised that some 
water exists in the liquid state as expressed by the unfrozen water 
content, Wie It has been shown that this phenomenon has relatively 
minor effects on ke and Ce in calculations for the rate of thaw. 


The remaining variable in the thaw solution that has yet to be 
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considered is the volumetric latent heat of the soil, L. This variable 
has a dominant role in thermal calculations and represents on a volume 
basis the amount of heat that must be supplied to a soil below 0°C to 
change ice into water. Problems arise in defining the value of L to 
be used. Firstly at the temperature, Ty» applicable in a problem the 
percentage Wh must be known. Secondly, the effect of a gradual change 
of state of ice to water over the temperature range Tg to 0°C must be 
considered. 

The first approach that can be taken is to calculate the amount 


of ice per unit volume once the appropriate Wa is known. That is, 
dyes Yq oll = Wa) c (2.515) 


where L' = 79.6 cal/gm the latent heat of water, 


Yq = dry density of soil (g/cm®). 


The rate of thaw can then be calculated from the Neumann for- 


mulation by using the total water content, w, to define Ke Cs fand Ce 


u 
and taking ke = O43 mcal/°C cms. This approach has the effect of 
"lumping" the total latent heat at O°C. That&is, ites assumed ‘that all 
the ice initially present in an unit volume of soil at Vy melts at 0°C 
rather than over the range Ty to 0°C. 


The five thermal variables Cy Ces k ke and L are all functions 


u° 
solely of water content when the soil is saturated. The complete solution 
to the Neumann problem presented in Fig. 2.2 requires the two temperatures 
Vy and Ty. and the parameters Cis Ky kes Ke Ke and L. As these six 
parameters are all functions of water content, it is convenient to plot 
them as shown in Figs. 2.7, 2.8 and 2.9. The conductivities ky and ke 


have already been given in Fig. 2.5. Therefore knowing Ties 'y and the 
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Fig. 2.7 Volumetric heat capacity of soils 
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Fig. 2.8 Thermal diffusivity of fine-grained soils 
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Fig. 2.9 Volumetric latent heats of soils 


200 


3] 


32 


water content of a saturated soil, all the thermal parameters necessary 
for an accurate solution to the Neumann problem in Fig. 2.2 may be 
obtained from the charts included. 

The effects of a gradual change of state from ice to water over 
the range Ty to OSG and the accompanying effects of temperature or phase 
dependence of frozen conductivity and specific heat are investigated 
in the next section, as a numerical solution is required to solve the 


heat transfer problem. 


2.3 The Temperature and Phase Dependence of Soil Thermal Properties. 


In the previous section it was assumed that all ice melted at 
the 0°C isotherm, and that the thawing soil could be divided into two 
discrete zones, each displaying a temperature-independent set of thermal 
properties. Since it has been well established that liquid water exists 
over a wide range of temperatures below Occ, it might well be argued that 
the thermal properties are dependent on temperature in the frozen soil, 
and that this should be accounted for more rigorously in the heat 
transfer formulation. Therefore, a more generalised equation for 
conductive heat transfer may be derived, and employed to establish the 
magnitude of any errors involved in the Neumann formulation described 
previously for soils. 

Consider the thermal equilibrium of a small layer 6x ata 


depth x. Taking heat flow as positive in the positive x-direction, 


(50), cl ObpstckOled stra et Lal (2.16) 


X 


where (SQ), = change in heat of the layer, 


(Q), 


heat flowing into the layer at x, 
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(Q), + 6x = heat flowing out of the layer at x 
(Q.), 


heat generated internally in 6x by water 
changing phase. 
In a time ét, the heat quantities in eq. (2.16) are replaced 


by their temperature dependent relationships giving 


s 90 90 
Deter AK Og Ot rm arath (Olina Obl oa. 
- L'w yy OW, 6x (2.17) 


where c_ (6) is the volumetric heat capacity of the soil-ice- 
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water mixture (cal/cm.~ ~C), 


k (6) is the thermal conductivity (cal/°C.cm.s), 


Le is the latent heat of fusion of ice (cal/gm water), 
Yq is the dry unit weight of the soil (g/cm?), 

W is the total water content (g/g), | 

OW is the fraction of w which changes phase in time st, 
Wh is the unfrozen moisture content defined as g.water/ 


g.(ice + water). 
Rearranging eq.(2.17) yields 


d Wh 
Fae eee deidt » (2.18) 
As Wh is considered to be solely a function of the temperature 
6, implying no hysteresis of the relationship, the following may be 


written 


d Wh d Wh 26 (2519) 


ay ts at MS 7 My | : 
Oy) 48 vee! son ir set 


nov ew ia xe at ertenaear 


baosfqsy 946 8) -p9 nt 2ote i p 
 privig sander ‘rina orud 6 


~s9f-froe ont to yi fosqso ‘doer AP amen a at 
| : eft") "sud eet etn 


« (2,m9.. .\I89) tiv Fsaubhoo emg 


' (197 Bw mo \ 69) sat To noraut 0 156 dha oat “ee 4 
Py! “mo \p) Froe ory +> Hig s¥mu eth ‘silt, “et in 

-(e\e) snetnoo wa el Fagor, ait at 

3 omba nt sesdq esprit]. rottws to natty ey ony h 


\"Stew.p 26 ant ab inate wiztom ri98 ie, 


34 


and bringing the last term on the right hand side of eq. (2.18) to the 


other side gives 


SL Ae et xan eee (2.20) 

The term in braces on the left hand side is often referred to as 
the apparent heat capacity of the soil. This term contains the straight- 
forward heat capacity of the $oi]-ice-water mixture, and the second part 
expresses the heat which is liberated or absorbed by the soil on freezing 
or melting. Hence the term 


dW 
ss 1 u 
CVC °P8 A) ca ONLY Sar or (2.21) 


is written as the apparent volumetric heat capacity of the soil. The 
quantity Cy (6) is written in terms of We from eq. (2.12), and therefore 


Cc. (9) becomes solely a function of Wye as follows 


dW 
U 


Using a method due to De Vries and described by Penner (1970), it 
may be shown that the temperature dependent conductivity may be well 


approximated by a linear dependence on Wh also, and hence 


k (@) = ke is Wh (k, - ke) (2223) 


where ky and Ke are the fully thawed and fully frozen conductivities 
respectively. 
Therefore all temperature dependent coefficients appearing 
in eq. (2.20) may be written in terms of Woe It remains now to establish 


the functional dependence of Wh on the temperature 86. 
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An exponential relationship between Wh and 6 is adopted for the 
purpose of examining the effect of temperature dependent thermal properties 
on the rate of thaw. It should be noted, however, that any continuous 
function might be used to relate Wo to 6 in this way. For example, 
Anderson et al., (1973) have expressed the unfrozen water content data 


versus temperature by the simple power equation 
Wa =a of 


_where a and 8 are constants which may be related to the specific surface 
area of the frozen soil. Dillon and Ancersland (1966) have also provided 
a relationship for the terminal unfrozen moisture content at temperatures 

below -5°C. However, for the purpose of observing the sensitivity of the 
| rate of thaw to the shape of the unfrozen moisture content curve, no 
regard need be given here to such correlations with more fundamental soil 
properties, and a simple empirical equation may be adopted. The unfrozen 
moisture is defined here in terms of three quantities P, Q and R which 
are capable of providing a wide variety of realistic relationships, 
that is 

w= (P +e? * Ry 100 (2.24) 


This equation and its temperature derivative may now be substituted 
into the expressions for apparent heat capacity and conductivity in eq. (2322) 
and (2.23). The governing differential equation (2.20) is written in 
finite difference form, and the details of the procedure for solving it 
are given in Appendix Al. The listing of the computer program and the 
accompanying function subprograms for conductivity and heat capacity 


are given in Appendix A2. 
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A set of unfrozen moisture curves covering a wide range are 
drawn in Fig. 2.10 by assigning values of P, Q and R to eq. (2.24). 

The values of P and R are constant for all the curves and are chosen so 
that Wh is unity at 0°C and 0.3 at the ground temperature Ty The 
value of Q is varied in order to produce the range of curves shown in 
Fig. 2.10. It is thought that this range of curves together with the 
Neumann solution (where all water melts at 0°C) must bound the 
relationships that are physically possible for real soils. 

The effect of these unfrozen moisture content curves on the 
movement of the O0°C isotherm is studied. Temperature dependent thermal 
conductivity and specific heat functions are established for the frozen 
zone, and the necessary parameters entered as data in the computer program 
given in Appendix A2. A simple example is solved where a step increase 
in surface temperature of 10°C is applied at the surface of a homogeneous 
soil exhibiting the three different Wy vs. temperature curves. The water 
content is 0.48, and therefore the conductivity of the thawed and fully 
frozen soils may be read from Fig. 2.5. 

The depth to the 0°C isotherm X is plotted for each case against 
the square root of time, and a straight line relationship is always obtained. 
Therefore, an a value may always be extracted from the solution. The 
results for X against /time are shown in Fig. 2.11, and the « values for 
the different Wh - temperature curves are compared with a Neumann solution 
in which the latent heat L is calculated by eq. (2.15) and "lumped" at 
0°C. Another Neumann solution in which Wh = 0 is also included for 
comparison and is shown by the dotted line in Fig. 2.11. 

It may be thought that the usual approach in dealing with the 


unfrozen moisture content relationship will underestimate the rate of 
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Fig. 2.10 Unfrozen moisture content curves for numerical analysis 
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Fig. 2.12 Effect of unfrozen moisture content curve on temperature profile 
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thaw. Jhis is true only if it is erroneously assumed that Wh =n §t 
can be seen that the rate of thaw obtained by the Neumann analysis actually 
overestimates by an insignificant amount the a value calculated when 
the melting range is taken into account. 

The agreement between a rigorously formulated model of heat 
conduction accounting for any reasonable phase dependence of conductivity 
and specific heat, and the Neumann solution presented on a single 
chart in Fig. 2.2 is remarkable, especially when the simplicity of the 
latter method is considered. There is no discernible difference between 
the a value calculated from the Neumann solution and that calculated 
from the first unfrozen moisture content curve in Fig. 2.10. The use 
of the second and third curves show very slight reduction in the a 
value. 

Although there is little deviation in the rate of movement of 
the 0°C isotherm, the temperature distributions in the frozen zone 
are certainly affected by the temperature dependent relationship. 
Temperature distributions for the same problems as solved in Fig. 2.11 
are presented in Fig. 2.12. 

The findings of this section might be usefully summarised at 
this point. Ina problem involving the thawing of a frozen soil,’ the 
dominant variables are the ground surface temperature, the thermal 
properties of the thawed soil, and the total quantity of water that 
changes state in a unit volume of soil. 

Two solutions are included for calculating an a value, the 
first being merely a graphical presentation of the Neumann solution 
in terms of three dimensionless variables. The second is a semi- 
empirical relationship which is accurate when the ground temperatures 


are close to zero. 
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The thermal properties of the frozen zone play an exceedingly 
minor role in the computation of thaw rates. It has been demonstrated 
that neither the absolute value nor the temperature dependence of the 
thermal conductivity of frozen ground is vital in calculating the 
depth or rate of thaw. For problems of geotechnical interest the use 
of the conductivity of ice equal to that of frozen soil will result in 
a sufficiently accurate solution. The volumetric heat capacity of 
frozen ground is easily estimated and small errors do not influence 
the thaw calculations significantly. 

It has been found that when soils exhibit a range of melting 
temperature the depth to the 0°C isotherm is stil] proportional to the 
square root of time for the boundary conditions of the Neumann problem. 
Moreover, it has been shown that there are no significant differences 
in the rate of thaw, a, calculated either by assuming that all ice involved 
in the change of state melts at aha ee or by accounting for the melting 
range of ice in soil. In fact the melting range retards slightly the 
propagation of a 0°c isotherm rather than accelerating it. A precise 
determination of the unfrozen water content - temperature curve is 
considered unnecessary in a thawing problem provided that the correct 
quantity of water changing state is input to the analysis. 

It can be concluded then that the solution obtained by lumping 
L as defined by eq. (2.14) is a satisfactory solution to rate of thaw 
problems. Ignoring temperature dependent latent heat effects introduces 
insignificant errors and, furthermore, is a mildly conservative 


assumption for soils. 


2.4 The Effect of Thaw Strain on the Neumann Solution 


As a thaw plane penetrates into a mass of frozen soil, considerable 
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movements of the surface boundary may occur relative to the fixed 
frozen soil zone for a variety of reasons, the most important of which 
are dealt with in this section. When an element of frozen soil thaws, 
there is a volumetric contraction of approximately nine percent in the 
water phase of the soil. If the soil is saturated, this effect will 
be responsible for an instantaneous volumetric strain of nine percent 


of the volume fraction of water in the soil, or 


Sa word Pel pe Mex (2525) 
where Say is the cumulative settlement due to the ice-water 
contraction, 


n is the porosity of the thawed soil, 

and =X is the depth of thaw. 

As the soil thaws, it may come under the influence of an externally 
applied stress, or perhaps just the influence of the self-weight of the 
soil itself. In either case, if pore water is present in excess of that 
normally in the soil pores under these stresses, the soil must consolidate 
by expelling excess pore fluids, causing a further strain on thawing. 

As will be shown in the next chapter, for many loading cases of practical 
interest, the degree of settlement, S/X, is a constant. It is considered 
useful then to combine the ice-water contraction and consolidation 
settlement in a relationship of the type 

Sue bak CARES) 
where 8B is the total strain experienced by the soil in passing 

from a frozen to a thawed consolidated state. 


This is done in order to assess the effect of thaw strains on the 


Neumann solution. 
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In general, the surface of application of the step increase in 
temperature moves relative to the frozen soil. This might also be 
the case in a situation where a hot pipeline was thawing the surrounding 
permafrost. If the bearing capacity of the soil were exceeded and the 
pipe were to sink into the soft thawed soil, the temperature source 
has moved closer to the thaw plane. 

These effects are most easily examined by considering the plane 
of application of the surface temperature to be the zero point of the 
Spacial co-ordinate x. This co-ordinate system places the frozen ground 
below in motion, with a velocity dS/dt which may be calculated from 


eq. (2.26) as 
eee ee 
dt dt (22h) 


The velocity v is negative if x iS measured downwards. 

A similar problem has been considered by Carslaw and Jaeger 
(1947) when accounting for the density change in the plane ice-water 
freezing problem. A solution is sought in which the position of the 


thaw line is given by 
X=avt (2328) 


and the temperature distribution in the thawed zone is given (as in the 


Neumann solution) by 


8, =A Grp eect on (2.29) 
2 


K 
u 


The equation of conduction of heat in the frozen zone moving with velocity 
v in the x-direction is derived for example by Carslaw and Jaeger (1947), and 
is 2 


0. 
ak Ce 2 a (2.30) 
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Employing eq. (2.27) and (2.28) in eq. (2.30) yields 


It may be verified that with the above value of X, eq. (2.31) is satisfied 


by 


@, = € erfc (~—— + a a ih ae (2.32) 


ay Ket av Ke 


The constants C and D may be determined from the boundary 
conditions, and the temperature distributions for 6, (x; t) and O5(x, t) 
may be shown to satisfy the governing equations and boundary conditions. 
The temperature distributions are not included here, but when substituted 
into the heat balance equation at the thaw interface, a may be determined 


as a root of the equation 


~My ae La (1 + Ba GI? 
SAGs eee al Ot) UI eg’ Ce 
k Ke me 
erf(—2—) u erfc [a (1 + B)/2/ Ke ] 
ear 
L/Yta 
av Ky cy (2.33) 


which is the same as the Neumann solution given by eq. (2.2), with the 
exception of the quantity (1 + B) appearing in the middle term. 

The thaw strain number B only appears in a quantity which is known 
to have a minimal effect on a, (see Fig. 2.2) when the ground temperatures 
are close to zero. It is difficult to see how the effect of B on the rate 
of thaw a might be expressed generally in graphical form, so it is 
proposed to perform two sample calculations of a at two different ground 


temperatures in the range of interest. We assune the following typical 
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set of properties. 


me eH pera 

w = 40% 

— 3’ cand 5°C 
Be emgyay: 


A thaw strain of 25% is certainly a large value which is not encountered often 
in natural soils, and the ground temperature of 56°C represents a lower 
bound to many engineering field problems. 

The calculations for a with a thaw strain of 25% and the specified 
ground temperatures were carried out using the revised equation (2.33), 
and the following results were obtained. 


FO, ae porsss a = 0.0328 cm/s!/? 


and for T, = -5 3 a = 0.0313 em/s!/@ 


The same calculations were carried out using the Neumann solution 
given in Fic. 2.2 which does not account for thaw strain, and these 


results were obtained. 


Sperone aso cy a = 0.0329 cm/s '/@ 


g 


and for T 
atid 


The deviation from the ordinary Neumann solution given in Fig. 2.2 


0.0314 cm/s //4 


It 

1 
oO 
e 

tt 


when a large thaw strain is accounted for is, therefore, approximately 
Lin ou0k or 02o7. 
It is concluded that in the grcund temperature range of interest, 
only negligible effects are observed when large thaw strains are 
accounted for theoretically. Consequently, provided the thaw depth, X, 
is measured from the point of application of the surface temperature, 
the effects of large strains in the thawing zone may be ignored completely. 


This assumption was implicit in an analysis carried out by Brown 
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and Johnston (1970) when calculating the depth of thaw under dikes 
constructed on permafrost at Kelsey, Manitoba. In this instance, the 
Strains in the thawing soil were known to be approximately 30%. The 
analysis assumed a linear temperature profile in the thawed soil and 
ignored the temperature distribution in the frozen soil. On comparing 
the depth of thaw (measured from the ground surface) with observed 


thaw depths, excellent agreement was obtained using this simple analysis. 


2.5 The Effect of Water Migration on the Neumann Solution. 


In the previous section, the effects of large strains in the thawed 
zone were examined and found to be negligible. Although a considerable 
quantity of water must pass through the thawed soil to produce large 
thaw strains, the possible transfer of heat by the melt water migrating 
under an excess pore pressure gradient has not been considered. This 
component of heat transfer, which is sometimes known as advection, has 
not been incorporated in an analysis for thawing soils. 

The equation of heat transfer in a thawed soil, including bulk 
motion of the pore fluid may be derived as follows. The derivation may 
be considered a coupling of the continuity of mass in the water phase 
and the equation of conductive heat transfer in a thawed saturated soil. 

The conductive heat flux at x is 

q = -Ky ox (2.34) 
where 6(x, t) is the temperature distribution in the thawed soil. 


The convective heat flux is 


Q=c, 6 v(x, t) (2535) 


where Co is the volumetric heat capacity of water, 


and v(x, t) is the velocity of the pore fluid defined on the basis 


of total cross-sectional area. 
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The heat transfer equation may now be written as: 


oe 2 eee ad 

at (Cy 8) Bae 9p) ian 8 (2236) 
or 

oC 
ele) OV ele) 3 ele, 
BDSG E aaa \ cerita a DED SE ae es 

Sta eS Cy 83x ~ Say Vogt ag Ky ox! (2.37) 

Since 

Ct va (0: 17) 4° ——w) 

W 

then 

oc Vv Cc 

SOULE eek: Sete iinet Sa te AS bape (2.38) 

ot yw w 9gt yw G. ot 
where e is the void ratio, 

w is the water content, 

G. is the specific gravity of soil solids, 

Yq is the dry density = Gy /O +e). 

Replacing Yq by its usual representation, we obtain from 
eq. (2.38) 

oC 

ue-¢ ] & 

ot. -.1te 2 (2.39) 

Now the equation for continuity of fluids ina soil may be 
written as 


ie teu Btu) if Uex (2.40) 


Pr) ru yas 
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and therefore from eq. (2.39) 


ee Ue L oe > OV 
5 CA itice at = Cw ox (2.41) 
° e ° ° oc 
Substituting this relationship in eq. (2.37), the terms in se and 
av 
dX cancel out, and we obtain 
Como Cie 90 
Beat. 5x (ky a CW v(x, t) ) X (2.42) 


and if Ky may be considered independent of position in the thawed soil, 


then 


rear ro Say t) = (2.43) 

Equation (2.43) may be used to predict the thaw depth for an 
equation of consolidation in thawing soils, and in turn, the consolidation 
equation provides the pore fluid velocity v(x, t). Therefore, the heat 
transfer and consolidation equations might be solved numerically for 
arbitrary conditions of thawing, loading and drainage conditions in a 
thawing soil. 

If an analytical solution to eq. (2.43) is desired for a step 


increase in temperature applied at x = 0, severe restrictions are placed 


on the form of the velocity function v(x, t). If a solution of the form 


X=avt (2.44) 


is sought, it appears to be possible only if 


dX 
v(x, t) Oc a (2.45) 


and it is doubtful if the case of a constant velocity might be solved. 
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Because of the obvious advantages in obtaining a solution in 
closed form, a velocity function of the type expressed by eq. (2.45) 
is adopted. It will be seen in succeeding chapters how these conditions 
for pore fluid velocity are achieved in many cases when a soil is thawing 
and expelling a constant quantity of excess pore fluid per unit depth. 
(See section 5.5 on the thawing of incompressible ice-rich soils). 

Letting the excess ice content, expressed as a strain, be denoted 

by E., the pore water velocity is given by the product of excess pore 


fluid content and the velocity of the thaw front, that is 


: dX 
vix, t) = - Es ge (2.46) 


where E. is the constant thaw strain due to the expulsion of 


excess pore water. 


Substituting eq. (2.46) into the governing eq. (2.43) provides 


00:5 v0 dX. 90 
X 
ai 
where..cekaca= op at and is dimensionless. 


A solution to the governing equation (2.47) in the thawed soil 


is sought in the form 


@=A+Berf (—~ ‘sonal eS (2.48) 
/ Kut av Ky 


which ensures that the governing eq. (2.47) is satisfied. On obtaining 


the constants A and B from the boundary conditions of the thawed zone, 
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the temperature distribution is given by 


(T. - Tp) ¢erf (—L_+k poral als “~— ) 


eV et ee OG 
Q=T - : u 
erf (—— + K See event (ih Seer) 
2/ Ky 2y Ky Oo Ky (2.49) 


where tT. and Te are the surface and freezing temperatures respectively. 


The rate of thaw parameter, a, is determined as a root of the 


transcendental equation 


The equation (2.50) has been derived ignoring the temperature 
distribution in the frozen ground. If Ty is not close to zero, then the 
middle term in eq. (2.33) may be included in eq. (2.50) to incorporate 
the effects of the temperature profile in the frozen zone. However, 
these effects have a minimal effect on a, and it is desired to isolate the 
effects of water migration in this study. Equation (2.50) is similar in 
form to the Neumann solution with UG = 0, which has been discussed earlier. 
One new dimensionless quantity Ky is introduced however, which indicates 
the effect of the advective heat transfer in the thawed zone. The 
normalised thaw rate af2v Kk, is shown plotted as a function of the Stefan 
number for two values of the ratio Se in Fig. 2.13. The possible range 
of values for K, must be between 0 and 0.5. As C /cy is usually around 
1.4, then K, = 0.5 implies that a thaw strain of 35% is occurring in the 
thawed soil due to the expulsion of excess pore fluids. As the value 


of a is affected so slightly in Fig. 2.13 by the possible ranges of Ky 
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Fig. 2.13 Effect of melt water migration on normalised thaw rate 
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it is thought necessary to show only the two curves for K, = 0, and 

K, = 0.5. For lower values of the Stefan number, the effect of 
advective flow must be considered to be negligible. At extremely high 
values of Ste greater than unity, the effect of a large K, value becomes 
noticable although still slight. 

It may be concluded that from the results of the parametric 
study given in Fig. 2.13, that the effects of edvective heat transfer 
in the thawed zone may be ignored in any practical analysis, without 
fear of incurring errors of greater than two or three percent. 

In addition, it may be noted that if the effects of water migration are 
ignored, ten any errors incurred result in a slightly conservative 
estimation of the rate of thaw. 

The other type of advective heat transfer which may be of 
concern is that due to lateral flow of water in a thawing slope or in 
the 'bulb' of thaw around a warm pipeline. If the flow is parallel to 
the slope of the ground, which is usually the case, then the flow of 
water in the thawed zone is parallel also to the temperature isotherms. 
Now,other than complex effects where the ground water enters and leaves 
the slope, ground water will not cross temperature isotherms, and 
therefore the heat transfer problem perpendicular to the ground 


surface is totally unaffected. 


2.6 Time Dependent Surface Temperature. 


It may often be necessary to investigate the velocity of a thaw 
front caused by surface temperature variations other than a step 
increase in temperature. In attempting to solve field problems within 


the framework of one-dimensional heat conduction, cases may arise in 
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which the prescribed surface temperature bears little relation to a 
step increase in temperature. 

The simplest method of introducing an arbitrary prescribed 
surface temperature into a method for calculating thaw depth in frozen 
soil is by the use of the Stefan solution. This solution assumes a 
linear temperature profile in the thawed zone, and ignores the 
temperature distribution in the frozen soil. The accuracy of the 
method has already been discussed in section 2.2 for the case of a 
step increase in temperature. If the surface temperature variation 
with time is denoted by the function T(t), then the depth of thaw 


may be shown to be 
2k eet dt 
X= /—A-? = — (2.51) 


where ie dt is the area under the surface temperature-time curve, 
and is known as the Thaw Index. 
This type of relationship is reviewed for example by Aldrich (1956). 
The Thaw Index is usually given in units of degree days. 
Equation (2.51) might be applied to determine the thaw depth under 
arbitrary specification of the surface temperature, with the same 
order of accuracy as the simple Stefan solution for a step temperature 
described in section 2.2. As an example, the depth of thaw under a 
sinusoidal surface temperature application is expressed as follows. 
Let the surface temperature vary as a half sine wave between time 0 
and time tes and denoting the maximum value by fae the expression 
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tec w-sin Ce (2.52) 


and : t 
t cena tix f Tt 
if uP dt = fi gas (1 - COS 075) (2.53) 


Placing eq. (2.53) in eq. (2.51) the depth of thaw is 


expressed as 


(2.54) 


As shown in Fig. 2.3, the Stefan solution for the step increase 
in surface temperature introduces significant over-estimations of the 
thaw rate a as the Stefan number increases. The Stefan number is 
defined as 3 Terk and the solution given by eq. (2.51) assumes that 
the "sensible heat" Cy ie is small in comparison with L. Therefore a 
steady state or linear type of temperature distribution is valid. 
However, if Cy ps becomes larger, then some curvature of the thawed 
temperature profile results, and it is useful to consider the addition 
of an extra term to increase the order of accuracy. 

An analysis which increases the order of accuracy in this way 
has been performed by Lock et al. (1969). Solutions are presented for 
the movement of the thaw interface due to the application of sinusoidal 
and power law expressions for the surface temperature. The analysis 
is subsequently extended by Lock (1971) to include arbitrary 


excursions of the surface temperature. The results for the sinusoidal 
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and power law temperature functions are summarised here, and a sample 
problem solved. 

If a sinusoidal surface temperature history is adopted, then 
the surface temperature for a thaw season of 180 days is written in 


a consistent notation as 


+ 
Seen. t 
= = Sin (} (2.54) 
max Cc 
where ee is the peak surface temperature, 
and t is a reference time of 180/7 days. 


The depth of thaw is given by 


X = <tnax Sc fc Jo sin (zt-)- 22 sin (xt-) Sin (¢-) ee 

An example of the use of this relationship is given later in 
section 5.2, and Fig. 5.5 presents the thaw interface history for a 
particular problem. If an equivalent step temperature is calculated as 
shown in Fig. 5.5, a thaw rate shown by the dotted line in Fig. 5.5 
results. The equivalent step temperature is defined as that constant 
surface temperature which would give the same thaw index over the 
thaw season as the surface temperature sine wave. The equivalent step 


may be shown from eq. (2.53) to be 
er ee Oe (22.50) 


The thaw interface velocities for the sinusoidal and equivalent step 
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temperature applications are also shown in Fig. 5.5, and are markedly 
different. The significance of this in geotechnical studies is 
indicated also in section 5.2. 

The solution for depth of thaw in a soil whose surface 
temperature varies according to a power law between 0°C and Tey over 


a time period of te is given by Lock et al. (1969) as 


ee ky Thnax Le "4 1/2 744m + 1)/2 
L n+ ] te, 


ne We auth pein elie nent 
BNSF (; + ) (=) (2.57) 


t 


where the surface temperature is described by 


i “ (: J 
A, Pepe Vege (2.58) 


and nis an arbitrary power. 


When n = 0, the solution for the step temperature increase is 
recovered. As pointed out by Lock et al. (1969), that if the surface 
temperature recrosses the freezing point of the material, then the 
final depth of thaw is the same as that predicted by the simpler 
Stefan solution given by eq. (2.54). Therefore, if one were interested 
solely in the maximum thaw depth at the end of a thaw season, the 


correct result is already contained in the simpler Stefan solution. 
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2./ Thawing in a Two Layer Profile. 


In many field problems of practical interest, the assumption 
that a deposit of frozen soil is uniform with depth may be entirely 
unrealistic due to the presence of a surficial soil layer having 
different thermal properties. Many undisturbed locations in Arctic 
regions exhibit an organic layer of varying thickness, overlying a 
mineral soil with a considerably lower moisture content. Examples 
of a surficial layer which has no latent heat associated with it might 
also be cited, such as a gravel layer, styrofoam sheet or pavement 
surface. In general, a layer of different thermal properties is 
introduced between the temperature application and the plane where 
melting occurs. 

No exact analytical solution exists for this problem, but 
bearing in mind the minor inaccuracies involved, the problem may be 
easily handled by the Stefan assumption of a linear temperature 
distribution in the thawed zone. Multi-layered soil profiles are 
handled in this manner by the U.S. Department of the Army (1966), 
but it is of interest to perform an analysis for a two-layer profile, as 
it represents the most common departure from homogeneity. 

Consider in Fig. 2.14 the case of a surficial layer of height 
H, thermal conductivity Ky» and latent heat Ly» overlying an infinite 
depth of soil with properties ko and Lo. Ignoring the temperature 
distribution in the frozen soil, it is then evident that the time to 
thaw the overlying layer completely due to the application of a step 


increase in surface temperature is given by 
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| L 
t = 
BemecksT (2.59) 
since : 2k,T. 
PANS ULS Ty 


from the Stefan solution for layer 1. 

When the depth of thaw X becomes greater than H, two layers 
of different thermal properties co-exist in the thawed zone. Assuming, 
for consistency with the Stefan solution that the temperature 
distribution in each layer is linear, a simple solution for the rate of 
thaw in the underlying layer may be derived. 


The heat balance relationships are: 


38) 985 
at x = H, ki 3x) 7 ke OK (2.61) 
30 
2 
and at Xe=oXUt), SND aga L, & (2.62) 


where on and °°2 are the temperature gradients in 


OX OX 
layers 1 and 2 respectively. 


This leads eventually to the depth of thaw 


Z 
ke Wott) k 
X = y ieee Cin ss dias FD, -(2- ) 4 (2.63) 
Ly ] 


and the interface temperature between the two layers T,(t) is given by 
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The melting temperature is always taken to be zero. 


A sample solution is included in Fig. 2.14 for the case of a 


60 


(2.64) 


peat layer overlying a fine-grained soil of 40% moisture content. The 


thaw rate through the lower material only is shown as a dotted line, 


and at large times the thaw rate becomes parallel to the thaw rate 


for the two layer system. The results of the simplified solution were 


checked rigorously by a finite difference solution, and the results 
found to be in good agreement. Linear temperature profiles were 
maintained in each layer throughout the thawing process. Seider and 
Churchil1 (1965) also note the maintenance of linear temperature 
distributions in a two-layer system of insulation and soil. 

The analysis may be applied also to calculating rates of thaw 
in the underlying soil when the surficial layer is dry, or without 
latent heat, as might be the case with materials such as gravel or 
styrofoam. In this case, as no latent heat "barrier" is present in 
the overlying layer, a linear profile of temperature is assumed to 
be established in the upper layer instantaneously. This assumption 
can be shown to be perfectly justified by considering the time 
factor Kt/H? for the upper layer. Experience with solutions to 
problems of the heat conduction type would indicate a time factor of 
0.5 to be adequate to establish an almost linear profile of 


temperature in the upper layer. The time to achieve this time 
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factor is 


t = 0.5 H¢/K (2.65) 


For a one foot layer of dry sand, for example, the time to establish 


a linear temperature profile would be of the order 
t = 0.5 x 900/2 x 107? = 2.2 x 10° seconds 


which is approximately three days. This time span is likely 
insignificant when considering the depth of thaw over a complete thaw 
Beason, or the life-time of an engineering structure on thawing 
ground. If large depths of gravel are under consideration, however, a 


complete numerical solution to the problem might be advisable. 


2.8 Thawing Around a Warm Pipe in Permafrost. 


Considerable problems are associated with the production of 
hot oil from vertical wells completed through deep deposits of 
permafrost. A thaw annulus is expected to grow to 50 or 60 feet 
within the first two decades of operation of such a vertical well 
(Palmer, 1972), (Koch, 1971). Just as dramatic are the thermal effects 
of a heated surface transportation pipeline on the surrounding perma- 
frost. 

Lachenbruch (1970) has presented estimates of the growth of 
the thaw 'bulb' around a horizontal hot pipeline in the Arctic. The 
stability of a hot pipeline buried in sediments that are under- 
consolidated when thawed is a cause of much concern. This is not 
because of the melting in itself, but because of the effects of the 


rapid rate of melting on the stability of the thawed foundation. 
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Although a fairly complex two-dimensional problem in heat transfer must 
be solved to obtain the rate of melting around the pipe, the most 
critical rate of thawing occurs vertically below the pipe axis. The 
relationship between thaw depth and time under the pipe may then be 
used in a one-dimensional analysis to determine the rate of 
consolidation of the soil under the centre-line of the pipe. 

The depth of thaw under the centre-line of a 48 inch hot pipe 
is plotted with time by Lachenbruch (1970) for two soil moisture contents 
and two ground temperatures. It is found that if the depth of thaw 
is plotted against (time)9°3, an almost perfect linear relationship 
is obtained in Fig. 2.15. It appears that the depth of thaw in this 


type of thawing configuration might be expressed as 


x = Bt? (2.66) 


where B_ is the thaw depth under the pipe centre-line 
after one year. 

If the hot pipe configuration envisaged by Lachenbruch were 
considered to be a one-dimensional thawing situation, the Neumann 
solution might be used to solve the problem. This is equivalent to 
applying a step temperature over the complete surface of the half-space, 
instead of the limited surface directly affected by the pipeline. This 
approach would obviously be expected to over-estimate the thaw depth 
under the pipe axis. However, for a short time at the initiation of 
thawing, the thaw depth should be well approximated by a one- 


dimensional analysis. As the depth of thaw becomes significantly 
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larger than the pipe diameter, two-dimensional effects become more 
important. Using the thermal properties and temperatures given by 
Lachenbruch (1970), an a value may be calculated from the Neumann 
solution, and the depth of thaw compared with the rigorous two- 
dimensional numerical calculations performed by Lachenbruch. The 
calculations for the one and two dimensional cases were carried out 
using a moisture content of 65% for 'clay', and 17 1/2% for 'silt'. 

The 'Arctic' location is assumed to have a ground temperature of -8.9 C, 
and -0.8°C for the ‘interior’. The temperature of the hot pipe is 

taken as 80°C. The one-dimensional a value represents the Neumann solution 
for these thermal properties (locations). The depth of thaw under a 
Surface temperature which is applied over an infinite Surface area 


is given as usual by 
) Make aha 


The corresponding B values for these locations are taken from the two- 
dimensional analysis given by Lachenbruch, where the thaw depth is 


calculated (as demonstrated in Fig. 2.15) from 


‘ ‘ 0.3 
At a time of one year after commencing operation, / t and t are 


both unity, and therefore the a and B values for the different 
locations given in Table 2.2 represent the depth of thaw in feet 
after one year. 

The depths of thaw using these one and two-dimensional 
representations are calculated also for a time of 0.5 years. Fora 


pipe with a diameter of four feet, Table 2.2 indicates that no 
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serious error is incurred by the one-dimensional analysis, providing 
the thaw depth does not exceed approximately 10 feet. When the thaw 
depth is at 14 feet, the one-dimensional analysis overestimates the 

correct thaw depth by 2 feet. At later times, the error in the one- 
dimensional solution becomes progressively larger as expected. 

The fact that the depth of thaw under the centre-line of a hot 
pipe can be well approximated by a one-dimensional analysis for the 
first half-year of operation is used when analysing the data from a 
test pipeline loop at Inuvik, N.W.T. The field results show that one- 
dimensional heat flow conditions are preserved at early times, and 
this proves extremely useful in the analysis of the initial (and 
usually the most critical) excess pore pressure conditions in the 
thawing pipeline foundation. These predicted and observed thaw rates 


are described in detail in section 6.3. 


TABLE 2.2 COMPARISON OF ONE AND TWO-DIMENSIONAL CALCULATIONS 
FOR THAW DEPTH UNDER A 48 INCH HOT PIPE 


1-Da 2-DB 1 - D Thaw 2 - D Thaw 
Location value value Depth at Depth at 


(ft/yr!/2) (Ft/yr?3) t= 1/2 year 1/2 year 


Arctic Clay 13.06 hes O.c3 9.0 
Interior Clay 13.92 13.0 9.84 10.0 
ArCENG 251 1t (eal ke) 16.3 14.94 13:0 


Intenior-siit 22.0) 17.4 16.03 14.0 
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CHAPTER III 


THE ONE-DIMENSIONAL CONSOLIDATION OF THAWING SOILS 


3.1 Introduction. 

In the previous chapter it is described how a solution to the 
problem of melting in a permafrost soil might be established. The 
movement of the 0°C isotherm with time may be determined from the heat 
transfer problem, and this isotherm is assumed to form the moving lower 
boundary to the region of interest in the associated problem of soil 
consolidation. The frozen soil does not enter into deliberations on 
soil consolidation, as it does not transmit large quantities of pore 
fluids, or participate in any significant deformations. 

On thawing, the soil becomes a compressible porous medium, 
characterised by finite values of compressibil?ty and permeability. 

The final effective stress in the soil may be calculated from a 
consideration of the total stress and the hydrostatic pore water pressure. 
The initial stresses in the thawed soil may be quantified by experiment, 
or as will be seen later, some rational assumptions may be made. 
Consequently the intermediate transient pore water stresses and 
effective stresses in the thawed soil might be described by a 
consolidation equation written in terms of excess pore water pressures. 

It is proposed to work for the present within the framework 


of the linear Terzaghi consolidation theory. This theory of course 
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involves several limiting and almost certainly inexact assumptions 
regarding the behaviour of the soil skeleton and the flow of water 
within the pore spaces of the soil. Frozen soil is a complex, non- 
homogeneous medium comprising at least three phases, and it is 
recognised that the subsequent consolidation of the thawed soil is 
probably a much more complicated process than the behavioural 
description that follows. However, the same might be said about the 
behaviour under load of natural unfrozen soil, yet it is encouraging 
to see the substantial power of simple theories such as the Terzaghi 
theory in accounting for a wide variety of phenomena in the field. 
Moreover, simplicity of a theory has the advantage that the results 
are often applicable to practice with greater ease. Subsequently, the 
limitations of the theoretical results are readily defined. Therefore, 
while cognizant of the limitations, in the following paragraphs the 
settlement of a thawing soil is explored using a well known theory of 
soil consolidation combined with a simple but important solution in 


heat transfer. 


3.2 Formulation of the Linear Theory of Thaw-consolidation. 

As illustrated in Fig. 3.1, a one-dimensional configuration 
is considered where a 'step' increase in temperature is imposed at the 
surface of a semi-infinite mass of frozen soil which is homogeneous and 
isotropic with respect to thermal and consolidation properties. The 
solution to the problem in heat conduction is presented in Chapter 2, 


and the movement of the thaw plane is given by 
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Fig. 3.1 One-dimensional thaw-consolidation 
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Other rates of thaw shall be considered in a later section. 
The lower boundary of the consolidating soil is now defined. 


The governing equation of consolidation is 


¢ Bu. du a6 
Voge ate at (72) 
where u(x,t) denotes the excess pore pressure, 


x denotes the depth measured from the ground surface, 
o denotes the total stress applied at the surface, 


and Cy denotes the coefficient of consolidation. 


If the applied stress does not vary with time, then eq. (3.2) 


becomes 
0<x< X(t) ; ful au, t>0 (3.3) 
; Cy ae ‘tee wie : 


The same Praesens inherent in the derivation of this 
equation when applied to other problems in soil mechanics will also apply 
to thawing soils. The linear void ratio-effective stress relationship 
for instance may introduce considerable errors if the soil is 
exceedingly ice rich. In this case, large volume changes must occur 
before significant increases in effective stress take place. 

The thawed soil is also assumed to be saturated, and 
therefore any solutions obtained will tend to overestimate the 


magnitude of the pore pressure if an air phase is present. 


Boundary Conditions 


We now consider the boundary conditions necessary to obtain 
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a solution to the consolidation problem. 
From eq. (3.1) it is obvious that initially, the thawed 


region does not exist, so the initial condition is 
X= OF thet sh¢6 (37a) 


The upper surface of the soil is taken to be free-draining, 
and so 


Wasa Us Or) ot 0 (3.5) 


At the freeze-thaw interface, a balance between the rate of 
liberation of excess pore fluids and their discharge from the thaw 
line must be established. 

If the soil is thawing very slowly and the discharge 
capability of the soil is large, then any pore fluids liberated will be 
expelled from the thaw line, and settlements will proceed concurrently 
with thawing. However, if a fast rate of thaw is maintained, and the 
discharge capacity of the thawed soil is not sufficient to expel all 
excess pore fluids, then excess pore pressures will be generated in 


the soil. 


As the soil is compressible, the boundary condition states 
that any flow from the thaw line is accommodated by a change in volume 
of the soil. Consistent with Darcy's Law, the rate of expulsion of 


water from the thaw line is 


dv _ _ Ak ou (x,t) (3.6) 
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where A is the cross-sectional area of an element, 


k is the permeability of the soil, 


ou (x, t) is the excess pore pressure gradient 
at the thaw interface, 


and in is the unit weight of water. 
In a small time increment At the thaw plane advances a 
distance AX. The volume of the soil element so thawed is 


V=AA X (327) 


The volume of water expelled from this soil element from 


Bqrno.0), 1S 


pk oll 3.8 
Y, 3X (Xt) At (3.8) 


AV 


and is equal to the change in volume of the element. 


Combining eq. (3.7) and (3.8) yields the volumetric strain in 


the soil element 


au 
avers, 655 HOS ©) (3.9) 
V x 

Yy dt 


For a compressible soil, the volumetric strain is related to 


the effective stress by a relation such as 


AV 


CO fiir I (3°10) 
V m,, Ao 


where Mm, denotes the coefficient of volume compressibility, 
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and 


Ao' is the change in effective stress corresponding 
to the volumetric strain. 


Combining (3.9) and (3.10) gives 


eho’ ou(X, t) 
Nor = V OX (3 


dt 


where the coefficient of consolidation is given by 


where 


and 


¢. = k (3; 


The total stress at x = X is 


o(X, th =P, +y Xx (3. 


ae is the surface applied stress, 


y is the bulk unit weight of the thawed soil. 


At the thaw line the pore pressure is the sum of the excess 


and hydrostatic components 


P(X, t) = u(X, t) + ¥,X (3. 
and so by definition the effective stress is 
OnMCXL theo Xy t) - Pi iXs t) = EA +! X= U(X, £) (3: 


where 


y' denotes the submerged unit weight of the soil. 


The effective stress increment at the thaw line is the 
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current effective stress less the initial effective stress, that is 


AOue= GX. tC) = oa (32 16) 


where on denotes the effective stress in the soil 
if no consolidation were permitted on 
thawing, and is called the 


restdual stress. 


Combining eq. (3.15) and eq. (3.16) we find 


on ee oe ye Ck t) (3217) 


and placing this expression (3.17) in eq. (3.11) gives the relevant 


boundary condition at the thaw line 


abexe= X(t): Pe - 05’ #4 X eee yh eee (3.18) 


There have been previous attempts to establish the boundary 
condition at the thaw line but they have been either incomplete 
(Tsytovitch et al. 1965) or incorrect (Zaretskii, 1968). Tne condition 
proposed by Zaretskii is similar to that developed here. However, 
in his derivation he takes a change in water content of the soil to 
be equal to the volumetric strain of an element and thereby introduces 


errors into his solution. 


3.3 Solution of the Linear Equation of Thaw-consolidation. 


The analytical solution to the linear Terzaghi consolidation 


equation (3.3) subject to the boundary conditions expressed by 
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equations (3.1), (3.5) and (3.18) is required. 

The essential feature to be contended with in problems of 
this type is the existence of a moving boundary. In the mathematical 
treatment it is necessary to satisfy a pore pressure condition at this 
moving boundary. For this reason, it might be anticipated that an exact 
solution for this type of problem could be difficult to obtain. For 
arbitrary rates of thaw, this is indeed the case, however the 
mathematical complexities are reduced considerably if thawing is 
considered to be proportional to the square root of time. 

There are two methods of obtaining an exact analytical solution 
(if one exists) to a problem of this nature. The first employs direct 
methods such as the separation of yariables technique or the Laplace 
transform, while the second employs indirect methods. An example 
of the latter is the semi-inverse method. A solution of a certain form 
is assumed and if, on evaluating the required constants from the 
boundary conditions, it is shown to satisfy all the necessary equations, 
it is the unique solution to the problem. An example of this method 
of attack is seen in the solution of the Neumann freezing problem 
in Carslaw and Jaeger (1947), p. 283-286. 

However, a more satisfactory method of obtaining a solution 
is by the more direct separation of variables technique. The 
boundary condition (3.18) contains the quantities (Fe - o,') swhich 
is a constant, and the term y'X, which is proportional to the square 
root of time. So it is necessary to derive a solution separately 
for each of these loading conditions, and then superimpose the results. 


Superposition of results is permissible here as both the equation 
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of consolidation and the boundary conditions are linear. 


Self-Weight Loading. 
Setting es - Gap) temporarily equal to zero, eq. (3.18) 


becomes 


Applying the transformation 


X 


“n= REC) 


to eq. (3.3), provides 


du _ 1 3u 

x X.0Z 
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The time derivative in the new (moving) co-ordinate system is 


transformed by 
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and so the time derivative becomes 


Q 
NI& 


C3221) 


Substituting these derivitives in the consolidation equation 


yields 
aby. Ay cal eae eee (3.22) 
ot X 57° yeah eau y 2 
The boundary conditions (3.5) and (3.19) become 
Zee Oise, US 0s. 20 >:0 (3223) 
oe tu 
a Es Goan ae vee 3 t>0 (3.24) 
X dt 
Separating the variables z and t requires a solution of the form 
uz, U)v=e F(z) 2 GC t) (3.25) 
It may be shown that this is only possible if 
G(t) =Vvt 
or ua= F(Z) TOE (3.26) 


Substituting (3.1) in equations (3.22) and (3.24) gives 
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By differentiating eq. (3.26), the following equalities are written 


as 
au. dF 
Ze OZ vt 
2 Z 
oe - of Tee 
OZ dz 
and oueome ih 
ee aye 


Substituting these derivatives into eq. (3.23), (3.27) and 


(3.28), we obtain in terms of F(z) 


2 2 
dF a dF 
ax, 5 OAS vA —-F?=0 (3.29) 
d2* 2 Cy dz 
Zz =20 F=0 (3.30) 
2 ¢ 
Eg ; ; a v dF 
sais ep atone le 340 Ge (3.31) 


Equations (3.29) and (3.31) are now reduced to ordinary 
differential equations, functions of the dependent variable z only. 


The complete solution to (3.29) in Gibson (1958) is 
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function term was erroneously omitted in Gibson (1958) ). 

The constants A and B are determined from eq. (3.30) and 
(3.31). From eq. (3.30), at z = 0, the terms inside the squared 
bracket are not zero, so to maintain F equal to zero requires 


the condition that 
A=0 (3.33) 


From eq. (3.32) and (3.33) dF/dz = B, and therefore from eq. (3.31) 
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Applied Loading. 


with depth, and temporarily setting y'X = 0,allows eq 


For the case of the loading hie - oe) which is constant 


written as 


toveq. (3.88) and (3.23). 


x = X(t); (Pe - ot) - U 


Transforming to z co-ordinates; 


z= 1 3 (P -~ oO 


A solution to the governing eq. (3.22) is sought subject 
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It will be found that the only value of G(t) 


possible to transform eq. (3.22) and (3.38) to ordinary differential 


equations is 


or 


‘thus 


and 


G(t) = 1 

u = F(z) 
ou . dé 
OZ dz 
‘Lal gt 
52° dz° 


(3.39) 


(3.40) 


tnefenod et dot dw Mea? lane — 
od O8 (BT-£) ips ewolls 0 = 


j2atenibto-o9 S oF eatmohedayl 


/ ; 
{SE} 0 4 i - ae c b> {* o <0 9) 
FR X Hy ‘ 


Joetdve Inpuoe 2? (S8.E) .pe sities aii OF oktutoe an 
(+a ¥o sutsyv yino ons Jeit bavot od) rh tw +1 AES. rds tng (8) 
fstinexotthb yientbyo of (86.8) lente ($5 &). «pS. 


Substitution of these derivatives in eq. (3.22) and (3.38) yields 
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Placing these values in eq. (3.43) and then in eq. (3.40) 


gives the solution 
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(3.41) 
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Again, the most obvious simplification is to set 


R= — in eq. (3.46) 
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Superimposing the results from eq. (3.36) and (3.47) finally gives 
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Dividing across by (Ps pica y'X) yields the normalised 


excess pore pressure 
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and z= UE) is the normalised depth to the thaw plane. 
Direct substitution of eq. (3.48) into the governing equation and 
boundary conditions shows that this solution does satisfy all 


requirements. 


settlements in a Thawing Soil. 


It is also of considerable interest to determine the average 
degree of consolidation in the thawed soil. For this class of moving 
boundary problem we define this as the ratio of the consolidation 
settlement that has occurred up to time t, to the total consolidation 
settlement that would occur if thawing were suddenly stopped at t. 
The settlement associated with the ice-water density change will be 
treated separately. 


The settlement ratio is 


Smax Xx 
0 {e, % es} dx 


where en er and e are the initial, final and current void 


(3.50) 


ratios respectively. 


As a linear void ratio/effective stress assumption has been 
adopted thus far, eq. (3.50) is re-written in terms of excess pore 


water pressures as 
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us 1 
St : A ae + Y'x aor u(x, t)} dx 
> (3050) 
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Rewriting the pore pressure equation (3.48) in terms of x, 
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Substituting the integral (3.54) in eq. (3.53) and again introducing 


Ws the expression for the settlement ratio becomes 
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(3.56) 
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So the settlement at any time t, Sys is readily calculated 


if the ratio St and § are known. 
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For the following two extreme cases of loading, the solution 


is Summarised. 
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3.4 Results of Solution. 

The complete solution to the problem of consolidation in 
thawing soils, subject to the assumptions made, is shown to be a 
function of the three independent variables z, We and R. W. is the 
ratio of the effective overburden pressure at the thaw line to the 
applied external loading. If an initial effective stress is present 
in the soil on thawing, it is subtracted from Por The thaw 


consolidation ratio R is a measure of the relative rates of generation 


and expulsion of excess pore fluids. 
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Normalized values of excess pore pressures have been computed 
from equations (3.47) and (3.37) for the two extreme values of W. and 
are plotted in Figures 3.2 and 3.3. Excess pore pressures for 
intermediate values of W. may be calculated from the more general 
eq. (3.49). It is of particular interest to note that for a soil 
consolidating only under its own weight the pore pressure distribution 
is linear for a given magnitude of the thaw-consolidation ratio R. 

The variation of the excess pore pressure gradient with R for this 
case is plotted in Fig. 3.4 in a form more convenient for use. The 
| degree of consolidation has been plotted against R for different 
values of We InFig, 345. 

For the cases where the soil is consolidating solely under an 
applied load, or solely under the influence of self-weight, the relation- 
ships describing the excess pore pressure distribution and degree of 
Beonsolidatian are all independent of time. A similar feature was found 
by Gibson (1958) when considering the pore pressures generated by 
deposition of material at a rate proportional to the square root of time. 
The dimensionless thaw consolidation ratio is a fundamental parameter 
which, although independent of time, has a role similar to the time 
factor in more conventional consolidation problems. 

As the excess pore pressure profile is time-independent, it 
is convenient to plot the excess pore pressure at the thaw line (z = 1) 
against R for the different We. ratios in Fig. 3.6, As the most 
critical pore pressure always occurs at the thaw front, this plot 
represents the most interesting facet of the solution. The curve 


for W. = 0 will prove particularly useful at a later stage when 
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PORE PRESSURE, 


0 Oi2 0.4 0.6 0.8 1.0 


Fig. 3.2 Excess pore pressures for the applied loading condition, 
W. = 0 


oy 7 \ ’ Tes. i me 2 its } rhe . E ree 7 
“outtibnos entbsof betfags sit aot eosueeend S09: 2290K3 ee) ree 
. any va OA | fae ee 


DEP = 
EPTH, oy 0H) 


PORE pressure, zt) 
LX 


0.4 0.6 


Fig. 3.3. Excess pore pressures for the self-weight loading 
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(dimensionless thaw consolidation ratio) 


Fig. 3.4 Variation of hydraulic gradient with R, bb = CO 
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Fig. 3.5 Variation of el owe with R 
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assessing experimentally the validity of the theory. 


3.5 Post-thawing conditions. 


If the movement of the thaw front were to cease abruptly 
at,say, time t = tes then a more conventional fixed boundary problem 
arises. For the case where the depth of thaw is proportional to the 
square root of time,. the position of the lower boundary Xe will be 
given by eq. (3.1). For this new problem, the governing differential 
equation (3.3) remains unchanged, however the initial values are given 
by the excess pore pressure distribution at the end of thawing. 

The boundary condition at the base is rewritten from eq. 


(3.18) to give 


ean er or ode (3.58) 
The velocity of the thaw front * is now zero, giving the impermeable 
condition. 


at x= Xe : — = (3.59) 


The formulation of the problem may be summarised as follows. 


Setting t, = 0 and Xe = L for notational simplicity, 


a 
su _ 3 
fo 10 Ce ee i Oe (3.60) 
: Vv 2 3 
OX 
Ce Oe: (es es x = 0 (ool) 
t > 0": oso ; x= L (3.62) 
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(P =o ') erf(R %) 
C= Oe >i = Bae e 5 ces are. Wi 7 oe Oe Ee S05) 
-R 1 + J 
erf(R) + © 2R- 
TR 


The analytical solution to this problem may be expressed 
in the form of an infinite series, however because of the complexity 
of the initial values in eq. (3.63), an extremely difficult integral 
emerges. It is possible, however, to obtain a solution for the self- 
weight (y' x) loading condition, as the initial values are linear, and 


it may be found from Carslaw and Jaeger (1947) to be 
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TT n= 0. :(2n, +41) (3.64) 
where Up = = 
if Tg ] 
ae} 
2R 
and is the excess pore pressure at the thaw line at completion of 
thawing. 
In dimensionless form, eq. (3.64) becomes 
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The excess pore pressure isochrones from eq. (3.65) are shown 
plotted in Fig. 3.7, and these isochrones are valid for all values 
of R during thaw. 

The post thaw settlement is defined as the settlement occurring 
after te divided by the total post thaw settlement taking place 
be tween te and time oo . The degree of post-thaw settlement for 
the post-thaw phase of consolidation may be obtained by integration 


of the infinite series (3.65) as follows: | 


2, 


1 
.T) n2_ .-MT 
I aeons \ SE Helis (3.66) 


The degree of post-thaw settlement for this case will be 


defined as 
ik ae dz 
ie 
S = ] ss aia CA a 
ie u(z,0) 4, (3567) 
Us 
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The initial values (u(z,0) are given by the triangular 


distribution 
u(z,0) = UpZ. 
DOK 
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Substituting the series expression (3.66) into (3.68) yields 


S21~2 > oye 
Ur ae (3.69) 


The degree of post-thaw settlement expressed by eq. (3.69) is plotted 
in Fig. 3.9 against the time factor. 

The post-thawing pore pressures and settlements for the applied 
loading ce, - o,') condition must necessarily be evaluated by a finite 
difference numerical method. Because the initial pore pressure 
distribution is always linear for the y' loading case, the degree of 
post-thaw settlement will always be independent of R. However, for 
the cea ~ a, loading case, the shape of the initial pore pressure 
profile is governed by R, and therefore the degree of post-thaw 
settlement will in this case be weakly dependent on R. A simple finite 
difference program was written to solve the Terzaghi consolidation 
equation subject to the initial values for the UR - Cees) loading 
case given by eq. (3.63). The excess pore pressure isochrones were 
then used to obtain the degree of post-thaw settlement for different 
R values during thaw. The excess pore pressures at the base of the 
thawed soil are shown plotted against time factor in Fig. 3.8. The 
degree of post-thaw settlement is plotted with the time factor in 
Fig. 3.9 for different R values. It is worth noting that little 
difference occurs in the settlement curves from R = 0.1] to R= 2. 


This is a broad range of values which bounds most R values likely to 


be encountered. 
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Fig. 3.8 Base pore pressures after completion of thawing, We = 0 
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Fig. 3.9 Degree of settlement after completion of thawing 
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The settlement and excess pore pressure curves for R = Go 
indicate that an initially uniform excess pore pressure existed 
at the end of thaw, and therefore the subsequent development of 
settlements and the dissipation of excess pore pressures will correspond 
exactly to the classical Terzaghi solution to this problem in soil 
consolidation. The most important settlement curve is that for 
R= 0.1. This is the same settlement curve which is obtained from the 
Classical Terzaghi solution for an initially linear profile of excess 
pore pressure, and is given by Terzaghi and Peck (1967) in their 
analysis for a hydraulic fill. This settlement curve is used extensively 
in the experimental verification of the theory as it is valid for all 
R less than unity (say), and for all loading conditions. It is 
interesting to note that the time factor for 50% post-thaw consolidation 
is 0.28, as distinct from 0.197 which is obtained from the classical 


Terzaghi solution for the oedometer boundary conditions. 


3.6 Discussion and Application of Solution. 


The degree of consolidation settlement and the excess pore 
pressures in a thawing soil have been shown to be primarily dependent 
on the R value. This dimensionless quantity relates a, which is an 
indication of the thaw rate, and therefore of the rate of liberation 
of excess pore fluids, to the coefficient of consolidation Cys which 
describes the ability of the soil to dissipate the excess pore fluids. 
When R is zero, no excess pore pressures are maintained, and 
settlements proceed concurrently with thawing. When the R value is 


large (i.e. unity), severe pore pressure conditions are maintained 
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at the thaw line, and the degree of consolidation settlement is 
considerably impeded. The R value may therefore be used as a 
convenient index as to the severity of pore pressure and settlement 
conditions in a thawing soil. When the excess pore pressure 
conditions are combined with a failure criterion of the Mohr-Coulomb 
type, the available shearing resistance of the thawing soil may 

also be calculated. 

Dwelling briefly on the two parameters contained in the R 
value, it is certainly safe to conclude that tne degree of uncertainty 
in obtaining a value for the thermal parameter a is much less than 
the uncertainty in estimating or measuring Cy The complete range of 
a values which one might imagine for a wide range of Arctic problems 


2 em/sec!/*, a total 


probably lies between 3 x 1072 and 9 x 107 
variation by a factor of three. Even accounting for possible errors 
in conductivity data and moisture content, one would not expect the 
error in the estimation of a for a particular problem to be greater 
than 20%. The coefficient of consolidation Cy» on the other hand, 
could quite easily vary from 107 cm?/sec for sandy soils to 107° 
cm@/sec for very fine-grained clay soils. Even though the Cy value 
enters the R value under the square root sign, its potential variation 
and the factors affecting it, make the degree of uncertainty 
associated with obtaining the correct Cy value much higher. Accepting 
that the Terzaghi equation correctly describes the pore pressure 
dissipation in the thawed soil, it is therefore encouraging to see 


that the solution of practical thaw-consolidation problems is resolved 


into a correct determination of the coefficient of consolidation. 
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The excess pore pressure conditions under the applied loading 
condition are in general slightly more severe than those obtained 
under the self-weight loading condition, but the difference is small 
enough that no generalisations should be made regarding the control 
of loading conditions at this early stage. It will be seen later, 
for instance, that the provision of an applied loading PS will in some 
cases reduce the normalised excess pore pressure in the thawing soil. 

The residual stress oe may be introduced as a constant or 
linearly increasing function into the preceeding analysis for linear 
thaw-consolidation without undue difficulty. 

Assume that the residual stress oe is peprecented as a linear 


function of depth by 


om (x) = D, + D, x (3.70) 


The boundary condition at the thaw line eq. (3.18) may then 


be rewritten as 


ou 
C eee 
(P, - D,) + (Y' - D,) x -u = ‘= (3071s 


t 


where the constant coefficient D, has been combined with the constant 
external stress Poe and the Do coefficient has been combined with y'. 
In the final solution, e is now replaced with (De - D,) and y' 

by (y' - Dy); and thus a residual stress function such as eq. (3.70) 
is immediately incorporated into the closed form solution for one- 


dimensional linear thaw-consolidation. 
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The theory formulated here provides a value of the settlement 
ratio Sak in the thawed soil. However, in order to calculate the 
actual transient settlements, a knowledge of the different components 
of total settlement and their relative rates is required. The total 
consolidation settlement na for an applied loading Re is proportional 
to the depth of thaw X, and therefore to the square root of time. 
ania for the self-weight loading case is proportional to e and 
therefore proceeds linearly with time. 

The other component of settlement arises from the volumetric 
contraction associated with the phase transformation from ice to 
water. This component of settlement occurs instanteously on thawing, 
and is obviously proportional to the depth of thaw and the square 
root of time, provided that the ice is uniformly distributed. 
Therefore in a situation involving both loading conditions, the total 
settlement is related to time and to the square root of time. For 
the oedometer conditions, however, where the self-weight component 
of the soil is ignored, the remaining components of total settlement 
are both related to the square root of time as follows. Writing 
the settlements for the oedometer loading conditions as the sum of 
the ice-water contraction and the consolidation settlements gives 

St 


S=J x +(e ee 
ma X 


where J is the strain associated with the ice-water change alone. 
The total consolidation settlement ee is written as 


max ; my, (P. fr 40 )x 
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Therefore the sum of the settlements at any time is 


( ee 
S = J+m PA Os X 
V a) 5 e 


For a uniform soil, the term inside the large bracket is a time- 
independent constant. Therefore the transient settlements taking 
place under oedometer stress conditions are related directly to the 

_ depth of thaw, and therefore to the square root of time. It will 

be seen in the following chapter whether these theoretical predictions 
are borne out for experiments carried out under oedometer loading 


conditions. 
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CHAPTER IV 
LABORATORY TESTING OF THAWING SOILS 


4.1 Thaw-Consolidation Tests on Remoulded Soils 

In order to assess the validity of some of the most important 
theoretical relationships derived in previous chapters, it was necessary 
to develop a special oedometer for laboratory testing of thawing soils. 
This oedometer was designed and built to test the theory of thaw- 
consolidation, and the results of the study are presented by Smith (1972), 
and subsequently by Morgenstern and Smith (1973). 

The oedometer was designed to satisfy the necessary stress and 
thermal boundary conditions during one-dimensional thaw-consolidation. 
The apparatus consisted basically of a rigid lucite thick-walled cylinder 
to maintain one-dimensional heat transfer, drainage and total stress 
conditions. Drainage is permitted at the top of the sample through 
a porous stone, and another porous stone at the base of the sample is 
connected to an electrical pore pressure transducer. A loading is 
applied vertically to the sample using a counterbalanced loading hanger, 
which bears on a plunger and load cap assembly. Temperatures at the 
top and bottom of the soil sample are controlled using thermo-electric 
elements and temperature control units built specially for the apparatus. 
Temperatures through the sample are measured using thermocouples and 
thermistors. The apparatus is described in detail by Smith (op cit), 


and a diagram of an improved apparatus appears in the next section. 
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Considerable care was taken in measuring the correct pore pressures 

at the base of the thawed soil sample. The pore pressures measured in 
tests performed on conventional unfrozen clay samples were found to be 
greatly affected by friction between the apparatus and soil sample, 
and by the 'hardness' of the pore pressure measuring system. The 
measuring system was made as incompressible as possible by keeping 

the volume of fluid in the measuring system as small as possible, and 
by carefully de-airing the fluid involved. Side friction on the sample 
was greatly reduced by lining the lucite barrel of the apparatus with 
a teflon sleeve, and a greased rubber membrane was placed between the 
soil and the teflon sleeve. In this way, almost ideal pore pressure 
responses were observed during standard oedometer tests on unfrozen 
clay soils. 

The soils were prepared by remoulding a saturated slurry at a 
high moisture content, and then by consolidation in the device to a 
desired stress. The remoulded soils were frozen in the apparatus to 
a uniform temperature. Upon the application of a sudden increase in 
surface temperature, the soil commenced thawing. Settlements, pore 
pressures and temperatures were recorded as thawing proceeded, and 
subsequently when thawing was completed. 

Using the curves for post-thaw settlement derived in Chapter 3, 
the value of the consolidation coefficient cy was calculated. When 
combined with the observed thaw parameter a, the thaw-consolidation 
ratio R was determined for each test. The settlement ratio, S4/S ray? 
and the observed pore pressure were plotted against R, and compared 
with the theoretical relationships derived from the thaw-consolidation 


theory. These results are reproduced in Fig. 4.1, and in this way a 
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Fig. 4.1 Observed and Predicted Pore Pressures and Settlements 
for Remoulded Soils (after Smith, 1972) 
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preliminary investigation of the theory of thaw-consolidation was 
accomplished. The results of the tests on remoulded soils, the 
measurement of the thermal and consolidations coefficients, and the 
analysis of the data within the context of the thaw-consolidation 
theory are described in a complete manner by Smith (op cit) and will 
not be dwelt upon here. The principal conclusions of this interesting 
study might be summarised in the following sentences. 

1) When a step increase in temperature is applied at the 
surface of a uniform soil sample, the depth of thaw penetration is 
observed to be proportional to the square root of time. 

2) The rate of thaw a may be estimated with sufficient 
accuracy for most engineering applications using well established 
relationships and thermal properties. 

3) Settlements during thawing under the application of a 
Surface step temperature are proportional to the square root of 
time, for oedometer stress conditions. 

4) The degree of settlement in the thawing soil follows the 
general trend predicted by the theory, but more precision in the data 
would be desirable. 

5) The linear theory of thaw-consolidation was more than 
satisfactory in predicting the excess pore pressures maintained 
during thawing. 

6) The versatility of the thaw-consolidation apparatus was 
demonstrated. From a single test, the parameters controlling the 
thaw and post-thaw settlements and pore pressures could be obtained, 
together with additional information on the rate of melting. 


7) The importance of the in-situ permeability for predicting 
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C, (and therefore R) was recognised when applying the theory to field 
problems. 

Following the extremely encouraging results of this study on 
remoulded soils, it was clear that further testing of undisturbed 
permafrost soils would be of considerable interest before applying the 
theory to field problems involving thawing soils. Although all tests 
carried out to assess the validity of the theory were performed under 
the applied loading condition in the oedometer, it is clear that any 
statements regarding the applicability of the theory apply equally 
well to the self-weight loading condition as the formulation of the 


consolidation problem is the same for each loading situation. 


4.2 The Design of an Improved Thaw-Consolidation Apparatus 

It was necessary to design and fabricate an apparatus capable 
of satisfying the stress and thermal boundary conditions involved in 
a thaw-consolidation test. The apparatus developed by Smith (1972) 
was shown to meet these requirements in a satisfactory manner when 
testing remoulded soils. However based on experience obtained from 
a large number of tests, and anticipating the extra complexities in 
testing undisturbed samples, Smith (op cit) provided the following 
recommendations for improvements to the original apparatus design. 

(a) The existing pore pressure measuring system was 
prone to plugging and leakage, and should be redesigned to allow 
flushing of the system. 

(b) The outer steel jacket in the original apparatus was 
not required for the low pressures involved in testing. 


(c) Provision should be made to measure deflection directly 
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from the load cap. 

(d) The sample should be mounted on a base plate, and a 
split oedometer barrel clamped around the sample. This innovation 
permits the placement of a frozen sample in the apparatus, which was 
previously impossible. 

These recommendations were carried out, and the improved design 
is given in Fig. 4.2. Except for the improvements listed above, the 
same features were incorporated in the oedometer used by Smith (op cit). 

Two thermo-electric modules powered by D.C. power supplies 
controlled the temperatures at the top and bottom surfaces of the soil 
sample. A lucite barrel lined with a teflon sleeve and a greased rubber 
membrane provided an almost frictionless boundary between the soil and 
the barrel of the oedometer. A lucite base plate was fitted with a 
porous stone and connected by a fine hole to the transducer and assembly 
of valves shown in Fig. 4.2. The system of valves allowed re-zeroing 
of the pore pressure transducer to atmospheric pressure without 
allowing the passage of any fluid in or out of the sample. An aluminium 
guide with a teflon bushing was provided for the plunger and load 
cap to prevent tilting of the load cap. A dial gauge accurate to 
0.001 in. mounted on the barrel of the oedometer allowed measurement 
of settlement directly from the loading plunger. The whole apparatus 
was mounted in a cold room, whose temperature was set at the required 
initial temperature of the frozen sample. The pore pressure measuring 
system was de-aired using a 50% solution of ethylene glycol to prevent 
freezing in the tubes or transducer. 


The method of measuring temperatures throughout the sample with 
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thermocouples and thermistors, and the method of data acquisition, was 
identical to that reported by Smith (op cit), and the reader is referred 


to this work for any further information on the instrumentation used. 


4.3 Thaw Consolidation Tests on Undisturbed Permafrost Samples 

The extensive laboratory testing carried out by Smith (1972) 
demonstrated very good agreement between theoretical predictions and 
the observed behaviour for remoulded soils. However, it is obviously 
desireable to extend this work by testing undisturbed natural samples 
of permafrost. It is extremely important to investigate the effects 
of the natural structure of permafrost on the predictions provided by 
the thaw-consolidation theory, and therefore determine the applicability 
of the theory to field problems. 

With the design improvements incorporated in the apparatus 
described in the previous section, it is now possible to place a frozen 
sample of the correct dimensions in the apparatus, and carry out a 
thaw-consolidation test. 

Samples of frozen core approximately 3 1/2 inches in diameter 
were obtained from Norman Wells, N.W.T. and from Noell Lake, Nortn-East 
of Inuvik, N.W.T. The core was inspected visually, and samples were 
selected for testing. The soil from Norman Wells was a sandy silt, 
with approximately 10% organic content. According to the NRC ice 
classification, the ice content varied from thirty percent visible 
ice (V. - 30%) at a depth of four feet, to about five percent visible 
ice (V. - 5%) from twelve to sixteen feet. Two small lengths of frozen 
core from the Noell Lake borehole were available from depths of five 


and fifteen feet. This hole was drilled close to the head scarp of a 
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large thaw subsidence which was occurring at the edge of a lake. The 
upper sample was a brown silty clay with twenty percent ice Ge - 20%), 
and the lower sample was a blue/grey silty clay with approximately 
twenty percent visible ice We - 20%). The ice in the upper sample 

was evenly distributed, and was present in the lower sample in the form 
of disoriented bands of 1/4 to 1/2 inch in thickness. Two photographs 
of this lower material are given in Fig. 4.3. A summary of the 
properties of these soils is given in Table 4.1. 

The samples were prepared in a large 'walk-in' cold room by first 
cutting a rough section of core to approximately 2 1/2 inches in length. 
The frozen soil was then clamped in a milling machine which was fitted 
with a specially constructed cutting tool. This cutting device was 
fabricated from a three inch long section of thin-walled steel tubing 
having an internal diameter of approximately 2.65 inches. The cutting 
edge of the cylinder was formed by bending the lip of the cylinder 
inwards and then machining the lip to a sharp edge. The internal 
diameter of the lip was machined to 2.5 inches exactly, and when 
mounted in the milling machine, the cylinder formed an effective device 
for cutting samples 2.5 inches in diameter for testing in the oedometer. 
With the cutting tool rotating at about 300 r.p.m., it was forced into 
the frozen core at a rate of 0.2 inches per minute. The cutting action 
produced a minute amount of melting at the cylindrical surface of the 
sample, but the effects were estimated to penetrate less than | m.m, 
into the sample. After the cutting tool had penetrated approximately 
two inches, it was withdrawn. The thawed "skin" around the sample was 
observed to refreeze almost instantly. The sample was removed from 


the milling machine, and the ends were trimmed square on a band saw. 
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NOELL LAKE 
CLAY 


DEPTH =15 feet 
Vs = 1.85 g /cm® 
w = 30-32% 

0= 0.2 kg/cm 


Fig. 4.3 Plan view of two samples of Noell Lake Clay 


. 


3#A) JU30 
Yaso 
igah@i= HT4a ao 
‘nat 2.1% 4% 
SRE-06 2 w 


‘mo\gv S.0 =o 


ian 


ysID stot Moot to 2otgnse owt Yo wabv wets Gb .p 
f 0 _. ate 


y 7 : : we 


5 o . oa r7G , an ” 
ie : j 
s* , % 1 a ‘ie th id Cpa A 
* ; ‘ae civ 4 


ee a > 7 - 


iad ORY 


7 
% 


114 


ALLS uoAag 


pole cul XS G9°z 22 o€ GS €2 poaynz Lysuoday 
Ae) 
a4e7 LL ON 
oe) coal “Ee GQ'°z €2 Lp 0S Ov paqunys Lpun 
Ke|) 
ADALY ULEZUNOW 
pol %2S'0 eal XL EL°2 = O¢ Ov a GG paznz LySsuod9y 
Ke) 
eoseqeuiy 
pol Gr 0 eal Mev GQ°Z 02 Ov vS Gt poqn} La suoday 
4LES 
SLLOM UeWUON 
oeUl - GeO! g-Ol X G2 GQ°Z 6L 62 L9 €2 pequnzs Lpun 
Wd wd 
S/ S/, (%) (%) aL eos 
AYLL Lqeauuad UOLZEPLLOSUOD 40 = Pw] 2 LWLT i eal ier 
LeoLdAy qualdijfao) LeoLdk] 9 IL4Seld DINDiAR Tals APL) % [LOS 


S3ILYAdOUd TOS 40 AYWWWAS L°p dav 


ak - 


¥ jf. 
1 136 


> 82 


a * JO 


$82 


tf 


20. 


= 
nd 


a ; 
7 a \ 
aye f 
: » 
ea 


3°49 


350 


WS 


The dimensions of the sample were carefully measured, and the sample 
weighed on a balance to 0.01 g. These data were used to calculate 
the frozen bulk density. The sample was quickly sealed in a plastic 
bag to prevent sublimation of the ice, and stored for later use. 
Trimmings and end cuttings from the sample preparation were used to 
obtain the initial moisture content of the frozen soil. Eight samples 
of permafrost were prepared in this manner, and the initial density and 
moisture content data are given in Table 4.2. Two reconstituted 
samples of Mountain River clay and Devon silt were also prepared to 
ascertain if the thaw-consolidation apparatus was functioning correctly. 
These samples were prepared by freezing de-aired slurries of the soil 
in lucite rings 2 1/2 inches in diameter. The samples were extruded 
from the lucite rings and trimmed for testing. The properties and 
initial data for these soils are also summarised in Tables 4.1 and 4.2. 
The sample was set up in the apparatus by first allowing the 
component parts of the apparatus to come to the temperature of the room 
at which the test was being carried out. As stated previously, the 
pore pressure measuring system and porous stone were saturated with 
ethylene glycol. A meniscus of glycol was maintained on top of the 
porous stone, and a wet filter paper was carefully placed on the 
glycol, and the frozen sample lowered gently onto the filter paper. 
The weight of the sample on the porous stone and base plate squeezed 
out all excess glycol, and a rubber membrane was placed around the 
sample. The membrane was sealed to the base plate by a rubber 0 ring 
seal. The two sections of the split barrel were clamped around the 
sample, and then the barrel was bolted to the base plate. The 


membrane was then held at the top of the barrel, and the load cap 
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assembly was pushed down the barrel to make contact with the soil 
sample. The overburden loading was then placed on the load hanger which 
was brought in contact with the load plunger. The vertical dial 
gauge was brought to bear on the load cap, and the initial reading 
was adjusted to equal the height of the sample. As the dial gauge 
had a 2 inch travel, this permitted the direct reading of sample 
height at any stage of the test. Any extraneous pore pressures in 
the measuring system were allowed to dissipate and the pore pressure 
measuring system was closed from the atmosphere. The assembled 
sample and apparatus were then left undisturbed for approximately four 
hours to allow the temperatures, height and pore pressure readings to 
become constant. 

The test procedure during the thaw-consolidation testing from 
this point onwards was the same as that described by Smith (1972). 
After applying a step increase in surface temperature to the sample, 
pore pressures, settlements and temperatures were recorded at one- 
minute intervals until thawing and settlement were complete. The 
results of the ten tests performed were plotted against the square 
root of time, and the results are shown in Figures Bl to B10 in 
Appendix B. The pertinent data are included on the graphs, together 
with sample data points. In order to calculate the settlement ratio 
S4/ Sia? ijt was necessary to calculate the height of the sample in 
the thawed, undrained state (Hi). This was done for each sample by 


assuming 100% saturation, and using the relations 
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and 


(4.2) 


where Orn is the thawed undrained void ratio calculated 
from eq. (4.'T), 
Y¢ is the frozen bulk density, 
and He and Hin are the frozen and thawed height of 


undrained soil. 


The excess pore pressures were normalised as in the theoretical 
treatment by dividing by ae ~ 5). P, was the applied loading, and 
the residual stress O5 in most cases was inferred to be zero due to the 
ice rich nature of the soil. In three cases, however, where the void 
ratio was low, it was measured for separate samples of the same material 
by the method described later in this chapter. The observed results 
for normalised excess pore pressures, settlement ratio and a values 
are presented in Table 4.2. The coefficient of consolidation was 
calculated from the post-thaw settlement data in the same fashion as 
that described in Section 4.1. When combined with the a value, the 
thaw consolidation ratio R was calculated. Figure 4.4 shows the 
normalised pore pressures plotted against the observed R value. The 
theoretical relation from Chapter 3 is also included as a solid line. 
The observed settlement ratios S4/S nex are plotted in Fig. 4.5 against 
the R value, together with the theoretical relationship. These results 
will be discussed in section 4.6, and are also given in Table 4.2. 


After complete consolidation had occurred, some additional loading 
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increments were applied in some cases to aid in defining the void ratio 
effective stress relationships, and obtain additional consolidation data. 

If a thaw strain parameter is defined based on total settlement 
from the frozen to the thawed consolidated condition as the total change 
in specimen height divided by the original frozen height, then the 
results of the eight tests on undisturbed permafrost samples may be 
compared with results by other researchers. Speer et al. (1973) 
produced a correlation between thaw strain and the frozen bulk density 
based on the results of many tests on fine-grained soils from different 
Arctic locations. This correlation is reproduced in Fig. 4.6 together 
with the results of the eight tests described here. The data points 
are observed to fall close to the statistical correlation given by 
Speer et al.(op cit) in most cases. 

The remaining study of importance is to determine the predictive 
power of simple analytical expressions for the thaw rate such as the 
Neumann solution, when combined with the simple expressions for the 
thermal properties described in Chapter 2. A value of a may be 
determined by calculating the heat capacities and conductivities from 
the data of Kersten (1949) and the known moisture content data. The 
latent heat term is calculated from eq. (2.15) based on the assumption 
that 10% moisture content (dry weight of soil) remains unfrozen at 
the initial temperature of -5°C. This assumption is based on observations 
from the published literature that the ranges of silts and silty 
clays tested most likely have unfrozen moisture contents at -5°C of 
from 5% to 15% of the dry weight of soil. Using these thermal 


properties together with the initial and surface temperatures of EHEC 
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and $219 Cs the graphical solution of the Neumann solution in Fig. 2.2 

is used to predict the a value. The observed value of a in the thaw 
consolidation test is then plotted against the corresponding predicted 
value of a in Fig. 4.7, for the eight tests on undisturbed samples. 

The dotted line at 45 degrees indicates the condition that the predicted 


a equals the observed a . 


4.4 The Residual Stress in Thawed Soil 

The boundary condition at the thaw line derived in Chapter 3 
equates any flow from the thaw line with a change in volume of the soil. 
In the linear theory, the volumetric strain of the soil is proportional 
to the change in effective stress. The change in effective stress at 


the thaw line Ao' is given by eq. (3.16) as 


Ao = o'(X,t) - o4 (3016) 


where o'(X,t) is the effective stress at the thaw line X; 
itpds a function. or, climes (ty, 
and om is the effective stress in the soil skeleton 


thawed under undrained conditions. 


The initial effective stress in soil thawed under undrained conditions 
will be referred to as the residual stress. Only departures from the 
residual stress will result in volume changes. 

The residual stress might conveniently be put equal to Zero, 
and this is a reasonable assumption if the soil is rich in ice, or has 
a high void ratio in the thawed, undrained state. However, it should 
be noted that in cases where the stress and thermal histories associated 


with the formation of a specimen of permafrost are such as to reduce 
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the void ratio of the soil, a significant residual stress may be 
present in the soil upon thawing. It has been shown in Chapter 3 
how the residual stress is introduced into the theoretical analysis 
for thawing soils, and a discussion of the physical aspects of 

the parameter may now follow. 

An appreciation of the physical meaning and importance of the 
residual stress within the context of the theory of thaw-consolidation 
is gained after consideration of the following experiment. A sample 
of unfrozen fine grained soil is prepared with a known stress history 
by consolidating a slurry to a known stress Pos The sample has 
dissipated all excess pore pressures, and therefore we is an effective 
stress (see Fig. 4.8). The sample is now frozen in the absence of 
drainage, and the average void ratio changes a small amount from point 
A to B due to the volume expansion of water to ice in the soil pores. 
If the sample is now allowed to thaw again without drainage, the void 
ratio returns of course to point A. However, a striking increase is 
noticed in the pore water pressure, and if the sample has a high 
enough void ratio, the pore pressure may approach or even equal the total 
stress on the sample. This implies that the effective stress is 
greatly reduced, and indeed may be zero for all practical purposes. 

If the sample is now allowed to drain freely, consolidation will occur 
and the sample regains equilibrium at a much reduced void ratio, 
point C on Fig. 4.8. 

Externally the soil in this sample has undergone a net decrease 
in volume represented by AC, under a constant external stress. During 
freezing, it is known that negative pore water pressures build up in 


fine grained soils. It may often be observed that small ice layers, 
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Fig. 4.8 Stress Path in a Closed System Freeze-Thaw Cycle (Schematic) 
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lenses, and other discrete segments of ice form upon freezing, even 
when freezing has taken place in the absence of free access to water. 
It must be concluded then that as the total quantity of water in the 
sample has remained unchanged, the remaining elements of soil are now 
overconsolidated with respect to the external total stress. That is, 
the areas of soil between discrete ice segments have experienced an 
effective stress greater than Fa There is an abundance of evidence 
for the existence of negative pore pressures during freezing (e.g., 
Williams (1966), Wissa and Martin (1968)). In cases where a significant 
clay fraction is present, it is well known that the negative pressures 
can be exceedingly high. 

On thawing, the overconsolidated elements of the soil possess 
an effective stress greater than ar However, a quantity of free water 
is available on a local scale from the thawed ice segments, and the 
soil swells almost instantaneously to absorb the excess water in the 
macro pores that result from the thawing of the segregated ice. The 
stress path taken by the soil elements is shown by the dotted line ADE 
in Fig. 4.8. If the soil is capable of absorbing all the free water, 
then the remaining effective stress is the residual stress. Alternatively 
if the soil swells to a zero effective stress condition, and free water 
is still available, then the residual stress is zero, and excess pore 
fluids remain in the soil. If free drainage is now permitted at the 
soil boundaries, the soil reconsolidates to the effective stress Poe 
and the reloading behaviour is typical of that of an overconsolidated 
unfrozen soil. The reloading path is represented by the line EC in 
Fig 4.0), 


The net strain from the frozen to the fully thawed consolidated 
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state (BC) is sometimes known as the thaw strain. In this simple 
experiment, the stress history of the soil is well defined, and it 
has been demonstrated how the thermal history is an extremely 
important factor in determining the residual stress in the thawed 
soil (5) and the subsequent stress-strain behaviour of the soil. 

The behaviour of permafrost is influenced by the stress and 
thermal histories, and the drainage conditions when the soil existed 
in a thawed state. When a sample of frozen soil is removed from the 
ground, the first measurement in the thawed state that may be made 
is the determination of the residual stress. This constitutes the 
starting point for estimations of the excess pore pressure and 
settlements in a thawed soil. 

If the soil is to remain essentially undrained for any 
Significant period of time for reasons of fast thaw rates, low 
coefficient of consolidation, or the lack of a free draining boundary, 
then the residual stress om will control the initial undrained shear 
strength of the soil mass. The relationship controlling the undrained 


strength of a purely frictional thawed soil is 


cy ‘ sect KC Ko) sin ¢ 
om he Hel2haly sina! 


where CC. is the undrained shear strength of the thawed soil, 
K. is the ratio between the lateral and vertical 
effective stresses under conditions of no lateral 
yield, 


A is the pore pressure parameter, 
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and ' is the effective angle of shearing resistance 


of the soil mass. 


In a soil which is simultaneously thawing and consolidating 
in a one-dimensional configuration, the excess pore pressures are 
controlled by the stress increment ie + y'X - of). If the combination 


of stress and thermal histories were such as to generate high values 


OF 6 
0 


» 1t is conceivable that O5 could be greater than ie + y'X), 
and negative pore pressures would result on thawing, accompanied by 
settlements in this instance smaller than those due to the ice-water 


transformation alone. 


4.5 The Measurement of Residual Stress 

In the following, an apparatus is described that will permit 
the measurement of the residual stress in a thawed soil. The 
restrictions that must be placed on the test configuration are that no 
water be allowed to enter or leave the sample during thawing, and 
that no lateral yielding be permitted so that one-dimensional 
conditions are ensured. 

The test might be carried out in either of two ways with a 
saturated soil. The first method is to set the total load o ona 
sample equal to some constant value (e.g., the effective overburden 
load that the sample might be subjected to in the field), and then 
after complete thawing the residual stress is calculated by measuring 
the pore water stress, and subtracting it from the total stress. A 
second procedure, which might be described as a ‘null’ method, can be 
undertaken by thawing the frozen soil and continuously adjusting the 


total stress o so that the pore water stress is always zero. In this 
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way, the residual stress on completion of thaw is equal to the total 
stress on the sample. 

The second of these methods is probably to be preferred, as 
it eliminates all pore pressure response effects in the soil and in 
the pore pressure measuring system. 

An apparatus of the oedometer type was designed and constructed 
to accept a 2.5 inch diameter sample of frozen soil. The principle 
feature is a split lucite barrel which was bolted to a lucite base 
(as shown in Fig. 4.9). A greased rubber membrane was placed around 
the sample with porous stones and filter papers at both ends. The 
membrane was then sealed to the base plate and the aluminum load cap 
with rubber '0' rings. The porous stones at each end of the sample 
were connected to a pore pressure transducer by means of the system 
of valves and tubing also shown in Fig. 4.9. This enabled the 
measurement of pore water pressures at the top and bottom of the sample 
at completion of thawing, and permitted the observation of subsequent 
consolidation behaviour of the thawed sample under additional load 
increments. The complete pore pressure measuring system was saturated 
with 50% ethylene glycol prior to testing to prevent damage to the 
valves and transducer in a cold environment. Temperature sensing 
equipment and temperature control devices are also present on the 
prototype, but these are not necessary for the measurement of the 
residual stress and the subsequent consolidation properties. 

The test was carried out by carefully machining the frozen 
sample to the required diameter, and trimming to a height of Leo 
1.5 inches. The sample was mounted in the apparatus as described 


above, and allowed to come to thermal equilibrium in a cold room at 
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at about -5°C. Great care was taken to ensure that the valves, tubes 
and porous stones were completely deaired with ethylene glycol. When 
the sample, membrane, and base plate were assembled, the two halves 
of the split barrel were clamped together around the sample, and 

were in turn bolted to the base. 

A large load of the order of 1 kg/cm® was placed on the frozen 
sample for a short period of time in order to ensure that the sample 
and load cap were seated. Any extraneous pore pressures in the 
measuring system were allowed to dissipate, and the valves were closed 
in preparation for the measurement of the pore pressure. An initial 
vertical dial gauge reading was taken, and thawing was commenced in 
an uncontrolled fashion, usually from the top downwards. It was usually 
adequate to thaw the sample by circulating warm water around the load 
cap, finally allowing the sample to come into thermal equilibrium 
with a room temperature of about +10°C. While the test was in progress, 
the pore pressure at the top of the sample was continuously monitored 
and the total load varied to keep the pore pressure zero. When 
thermal equilibrium was regained at the positive temperature, the total 
stress on the sample was then equal to the residual stress in the 
thawed soil. 

Additional load increments were then placed on the sample to 
determine the permeability, compressibility, and consolidation 
coefficients in the usual manner. 

On completion of testing, the sample was removed, its water 
content determined, and the void ratio at the thawed undrained 


condition an was Calculated. 


Measurements of the residual stress were made on reconstituted 
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samples of Athabasca silty clay from Northern Alberta, and on samples 
of a blue silty clay from the Mountain River, N.W.T. A summary of 
the properties of these materials is given in Table 4.1. The soils 
were prepared as slurries, and deaired. In some cases the slurries 
were consolidated to a convenient void ratio and then frozen in the 
apparatus. In other cases, frozen samples were prepared from the 
slurry, and after machining to the correct dimensions, were placed 
in the apparatus. 

In either case, the sample was thawed, and the residual stress 
measured as described above. At completion Of thawing, the void 
ratio of the sample was compared with the void ratio prior to freezing. 
For a successful test, these void ratios should be the same, demonstrating 
that no drainage occurred through the freeze-thaw cycle. Additional 
load increments were placed on the thawed soil and the consolidation 
properties defined. The sample could then be refrozen without drainage, 
and the residual stress measured for the lower void ratio. This procedure 
could be carried out several times until the relationship between the 
thawed void ratio e and om was defined. A typical set of test results 
for Athabasca clay, series A4, is shown in Fig. 4.10. It is worth 
noting that any relationship between e and oe does not represent a 
stress path taken by the soil, but rather the locus of a series of 
points taken from a set of stress paths such as those indicated in 
Fig. 4.10: 

In the test series A4, the soil was remoulded and consolidated 
to C.a2 kg/cm. The sample was then frozen for the first time, and the 


residual stress was found to be immeasurably small. The soil was 
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consolidated to a higher stress and frozen. Upon thawing the residual 
Stress was measured, and the process was repeated several times after 
consolidating to a higher load each time. To determine the effect of 
a different stress history on the residual stress, another remoulded 
sample was consolidated to 1.0 kg/cm? in test series A6, shown in 

Fig. 4.10 also. It was then frozen for the first time, and the 
residual stress measured in the thawed soil. Subsequent measurements 
Of o were undertaken by reconsolidating the sample to the same stress 
Ore ls0 kg/cm. 

As freeze-thaw cycling continued in the absence of drainage 
during freezing, the residual stress increased steadily until it 
approached the overburden stress. The consolidation which occurred 
on reloading steadily reduced as the cycling proceeded. This behaviour 
is to be expected, as it is sntoncetatere that continued cycling would 
continue to cause steady decreasesin the void ratio ad infinitun. 
Furthermore, this test series demonstrates that if the effective 
overburden stress is close to the residual stress in the thawed soil, 
then subsequent settlements due to consolidation after thaw will be 
small. The combined results for tests A4 and A6 are plotted in Fig. 4.11 
as thawed void ratio against the logarithm of the residual stress ome 
Although entirely different stress paths were taken by the two soil 
samples, approximately the same relationship is found between e and ome 

A strong linear correlation exists between the void ratio 
of the thawed soil and the logarithm of the residual stress. This 
relationship appears to be sensibly independent of the previous stress 


and thermal history effects, at least for the limited studies undertaken 


so far. 
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Two remoulded samples of Mountain River clay were subjected 
to different sequences of consolidation and freezing. The stress 
paths taken by the two test samples MR1 and MR2 are shown in Fig@4.122 
A summary of the residual stress data is shown in Fig. 4.13. Again 
the uniqueness of the void ratio-log residual stress line is demonstrated. 
A summary of other typical consolidation properties obtained during these 
tests is given in Table 4.1. 

Undisturbed samples of Norman Wells silt were from different 
depths, cut to the required dimensions and a profile of void ratio with 
depth was obtained from moisture content determinations. The frozen 
bulk density of each sample was measured before testing. 

A sample was placed in the apparatus and allowed to come to 
a uniform temperature of =5°C. | The sample was then thawed as described 
previously, and the residual stress was measured. The results of tests 
on several samples are given in Fig. 4.14. To obtain some additional data, 
further consolidation to a lower void ratio was allowed under additional 
loading, and the sample refrozen prior to measuring the residual stress. 
Although the samples were no longer undisturbed, the data are included 
in Fig. 4.14. The depth, void ratio, frozen bulk density, residual 
stress, and coefficient of consolidation are summarised in Table jose 

Again a linear relationship between void ratio and the logarithm 
of om is found. Fig. 4.14 also shows the virgin compression curve for 
a remoulded sample of the material not subjected to any freezing. 
Generally, the increased scatter in the residual stress data for the 
samples at higher stresses simply reflects the increased difficulties 


jn measuring pore pressures as the soil becomes less compressible. 
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It appears that for the site in question the residual stress 
increases approximately linearly with depth. The residual stress is 
shown plotted with depth in Fig. 4.15. The variability of the test data 
between 14 and 16 feet reflects the local changes in void ratio at 
that depth. 

It is clear that considerably more testing would be required 
to establish the relationship accurately for a relatively incompressible 
material such as this. However, it is encouraging to observe that a 
material which possesses five to ten per cent visible segregated ice 
(see Table 4.3) exhibits a significantly large residual stress on 
thawing. 

In addition to the measurement of the residual stress, the 
coefficient of consolidation was determined for the thawed samples of 
undisturbed permafrost. This was done by placing a load equal to the 
effective overburden weight on the sample, and allowing consolidation 
to occur from the initial effective stress (co) to the final effective 
stress (Pa + y'x). In this way, the coefficients Cys Mm and k were 
determined accurately over the correct stress range. The results for 
the undisturbed samples of Norman Wells silt are shown as a graph of 
the coefficient of consolidation with depth in Fig. 4.16. It is seen 
that below the permafrost table the coefficient of consolidation is 
relatively constant at 2.5 x 10°cm?/s. This value is in all 


likelihood large enough to ensure a stable foundation at this site if 


thawing were to occur at any reasonable rate. 


4.6 Interpretation of Laboratory Test Results 
The laboratory testing carried out by Smith (1972) demonstrated 
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Fig. 4.15 Variation of Residual Stresses with Depth 
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that the rate of thawing through remoulded samples of frozen soil can 
be predicted with reasonable accuracy. It was also shown that the 
excess pore pressures and settlements in a thawing soil are dependent 
on the R value, and agree well with the predictions made by the linear 
theory of thaw-consolidation. 

Using the modified design of the thaw-consolidation apparatus, 
thaw-consolidation tests were successfully carried out on several 
undisturbed samples of two Arctic soils, and two other reconstituted 
soils. The agreement between the observed rate of melting and the a 
value determined from the Neumann solution is remarkable, with 
deviations usually less than 5-10%. Bearing in mind that no special 
consideration was given to the unfrozen moisture content relationships 
for the soils, and that the conductivity data of Kersten (1949) were 
used to predict the a value, the predictive power of this simple 
procedure is excellent and extremely encouraging. The two parameters 
a and Cy enter into the R value, which is known to control the 
transient behaviour of pore pressures and settlements in thawing soils, 
and if a can be determined in advance to the degree of accuracy shown 
in Fig. 4.7, then attention may be diverted from heat transfer problem 
and directed towards resolving the consolidation problem and obtaining 
the correct value of the consolidation coefficient Cy to use in analysis. 

Using the Cy value determined in the post-thaw phase of 
consolidation and the measured ao value, an R value is determined for 
each test. In Fig. 4.5 the theoretical and observed settlement ratios 
are plotted with R, and the results are generally in good agreement. 


The two test results showing very low settlement ratios are most likely 
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due to non-homogeneity of the ice features in the sample. The scale 
of the ice features in these samples (I2 and 13) were large when 
compared with the sample size, and in all probability, large ice 
layers were encountered by the thaw plane near the base of the 
samples. In Figs. B5 and B6, the settlement against root time curves 
are non-linear, demonstrating the probable non-homogeneity of the 
samples. In general, however, the ability of the theory in predicting 
the settlement ratio is reasonable. 

The normalised excess pore pressures when interpreted in the 
context of the linear theory demonstrate excellent agreement between 
prediction and observation. Considering the extremely transient 
nature of the heat and moisture flow in the soil, and the general 
difficulty with pore pressure measurement in any configuration, these 
results are felt to constitute a successful verification of not only 
the theory of thaw-consolidation, but also of the functioning of the 
laboratory equipment. In four of the eight tests on undisturbed 
material, the thaw strain in the soil was less than twenty per cent. 
In the remaining four tests, very much larger strains occurred. 

This large strain behaviour is best predicted in advance by using the 
frozen bulk density of the sample as in Fig. 4.6, and by essentially 
disregarding the ice features in the sample. This observation is 
based on the results from samples I-2 and I-3, where the visible ice 
amounted to 20% by volume. However, the thaw strain under an 
approximate overburden load was only of the order of 12%. This 
strain could have been predicted by use of the frozen bulk density 


correlation given in Fig. 4.6, whereas an estimation of strain based 
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on visible ice would have overestimated the strains for these samples. 

It appears from the results given in Fig. 4.4 that the excess 
pore pressures in the samples exhibiting large strains were 
considerably higher than those predicted by the linear theory of thaw- 
consolidation. A better correspondence between theory and observation 
is obtained by incorperating a non-linear void ratio-effective stress 
function to describe the high degree of non-linearity of the constitutive 
relation for the soil when large strains occur. This analysis has 
been carried out, and is described fully in section 5.3. The method 
of plotting the void ratio against effective stress, and the use of 
this data to determine oF is described also in section 5.3, and was 
used here in the interpretation of the four tests exhibiting large 
strains. The value of om was combined with the applied stress EA 
to obtain the stress increment, and then the normalised pore pressure 
curves in Fig. 5.7 were used to predict the excess pore pressures in 
the thawing soil. The effect of applying this non-linear e - o' 
function in interpretation of the excess pore pressures is shown in 
Fig. 4.17. The pore pressure when interpreted in the context of the 
linear theory is plotted, and then the resulting change in predicted 
pore pressure is shown when the non-linear theory is used to interpret 
the results. This extension to the linear theory is shown to improve 
the predictive power of the linear thaw-consolidation mode] 
considerably, and the predicted and observed pore pressures no longer 
differ greatly. 

Finally, the frozen bulk density appears to corre late 
reasonably with the thaw strain, and the relationship proposed by 


Speer et al. seems to fit the observed soil strains in a satisfactory 
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manner. Although the correlation is entirely empirical in origin, 
it appears to be capable of estimating the total settlement 
approximately, and is therefore useful in a preliminary analysis. 

The existence of a residual stress in permafrost thawed under 
undrained conditions has far-reaching implications. The higher the 
residual stress in the thawed soil, the smaller will be the subsequent 
consolidation settlement, the lower will be the pore pressures generated 
during thaw and the higher will be the undrained shear strength. If 
the depositional sequence and the thermal history associated with the 
formation of a body of permafrost are such that the void ratio decreases 
significantly with depth, substantial values of the residual stress 
would be anticipated. Hence, the thawed soil would have a finite shear 
strength and problems such as the stability of the thaw bulb around a 
buried warm oil pipeline become less acute. The residual stress wil] 
generate a finite frictional resistance and the thawed soil will not 
behave like a viscous fluid as suggested in the prognostications of 
Lachenbruch (1970). 

The influence of the residual stress is also of concern in the 
consideration of the effects of thawing permafrost on the behaviour 
of an oil well. The arching of the soil about the well will affect the 
stresses transferred to the casing which tend to make it buckle. It is 
of interest to note that Palmer (1972) has drawn attention to the control 
of arching by the effective stresses existing in the thawing soil during 
thaw and subsequent consolidation. He refers to the “extent of initial 
consolidation" which is related algebraically to the residual stress 


and observes that it is a property of the site rather than of the soil 
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type. He notes further that it can only be determined by thaw 
consolidation tests on undisturbed cores. 

In many cases of foundation design it may be prudent to remove 
the top few feet of highly compressible ice-rich material in order to 
limit thaw-induced deformations of the foundation. However, unless 
the residual stress can be shown to have a significant value compared 
with the effective overburden stress, serious deformations may still 
occur as the thaw plane penetrates deeper into the underlying ground. 
On the other hand, if a stratum with a high residual stress is 
accessible without excessive excavation, it may be possible to found 
flexible structures directly on it without concern for subsequent thaw. 
For example, oil tanks might be buried in the permafrost rather than 
placed on the extremely thick pads of gravel currently used... 

The existence of the residual stress enables one to develop 
Support conditions in thawing ground with greater confidence in their 
performance. This study of the residual stress has not only improved 
the ability to predict excess pore pressures and deformations in a 
thawing soil but it has also provided further insight into some of the 
processes involved in the formation of permafrost. Much more 
experimental work on natural permafrost samples is necessary to gair 
further knowledge concerning the distribution of residual stresses 


in the field and their effect on the behaviour of thawing soils. 
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CHAPTER V 


EXTENSIONS TO THE THEORY OF CONSOLIDATION 
FOR THAWING SOILS 


5.1 Introduction 

In Chapter 3 a theory of consolidation for thawing soils is 
developed. Although the simplicity of the treatment allows the 
expression of the solution in an analytical form, and consequently 
permits examination of the variables controlling the excess pore pressures 
and settlements, it is recognised that important deviations from the 
assumptions made for the theory given in Chapter 3 may occur. In the 
following section several of the more important and obvious limitations 
of the simplified theory are discussed, and extensions are provided to 
increase the versatility of the theory. The more comprehensive the 
extension of course, the more complex will be the resulting solution. 
At the present stage of knowledge, it is considered more important to 
formulate the simple physical statements associated with each extension 
in turn, and observe the effects of the most important parameters, 
rather than attempting to establish some general, all-purpose numerical 
procedure that would account for many extensions simultaneously. When 
ijt has been decided which effects are important, and which are likely 
to be encountered in practice, then the implementation of such a general 


procedure might be suggested. 


Each of the theoretical cases considered in this chapter is 
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designed to extend the range of practical problems that may be 


embraced by the theory set out in Chapter 3. 


5.2 Arbitrary Movement of the Thaw Interface 


One of the most serious limitations of the linear theory 
described above is the assumption that the thaw interface moves 
proportional to / t. This is a valid relationship for the important 
case where a step increase in temperature is applied at the surface of 
a homogeneous frozen mass of soil. But for reasons of non-homogeneity, 
varying surface temperatures or two dimensional effects, a more general 
relationship between thaw depth and time may emerge from the associated 
problem in heat transfer. In section 6.2 such a relationship is 
extracted from a solution for heat conduction around a hot pipeline in 
permafrost. 

An analytical solution in closed form is available for the 
case of thaw depth proportional to Yt, but it appears unlikely that 
such exact solutions might be obtained for arbitrary movements of the 
thaw front. However, numerical methods may be employed to obtain a 
solution for the consolidation problem to almost any required degree 
of accuracy. 


Let the movement of the thaw front be given by 
X(t) = Bt” (5.1) 
and its velocity by 


dX 


en pro! (5.2) 
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This power law permits expression of a variety of realistic 


relationships between X and t, with B and n specified as constants. 


Rewriting the boundary equation (3.18) at the thaw interface 


using equations (5.1) and (5.2) provides the set of equations 


teen Os u= 0; x= 0 C503) 
; OU s 9¢y 
tat): aber Soran O, Sa Xacee dil te) (5.4) 
OX 
au 
t 1p4n 
t> 0; (Py - of) + y'Bt?-u= “8X5 x = x(t) (5.5) 
n Bt 
n 


These equations may be written in finite difference form, using 
a scheme proposed by Crank and Nicholson (1947), and the moving boundary 
is accounted for by either of the methods described by Murray and Landis 
(1959). The advantage of the Crank-Nicholson scheme is that it remains 
unconditionally stable for all values of the time step At, and the 
accuracy is much improved over other somewhat simpler methods. 

The simultaneous equations so produced are tridiagonal when 
written in matrix form, and a simple algorithm derived from the 
Gaussian Elimination technique for simultaneous equations provides a 
fast efficient solution. The details of the finite difference scheme 
and the Gaussian Elimination algorithm are summarised in Appendix A.3. 

Excess pore pressures u(x,t), were evaluated at selected 
times by a program written for a digital computer. The degree of 


consolidation in the thawed soil was evaluated as in the linear theory 
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from the expression 


2 ee (5.6) 


The integral term in equation (5.6) may be evaluated using a simple 
numerical integration technique such as the trapezoidal rule. A 
listing of the computer program is presented in Appendix A.4. 

Values of the normalized excess pore pressure at the thaw 
interface, and the degree of consolidation are plotted against a time 
factor for the two extreme loading conditions y' = 0, and (Ps - on) = 0 
in Figs.-5.1, 552;-5-3-and-5:4.--“1t-was-found- that-for-d=givenvalue 
of the power n, a unique relationship emerged between excess pore 
pressure, or settlement, and the time factor. The time factor T is 


defined as 


cyt 


i Hegre), 


~ x(t)? 
Values of n from zero to unity were chosen to cover the full 

range of physically realistic problems. A value of n equal to zero 

corresponds to a conventional fixed boundary consolidation problem, 

with the depth of consolidating soil set permanently equal to B, and 

an initial excess pore pressure equal to ues - oo) + y'X. These results 

for n = 0 are plotted for comparison purposes, to form a lower bound 

to the family of curves. 


The constant velocity thaw front (n = 1) is thought to be the 


upper bound for this set of curves, as this thaw rate is only achieved 
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Fig. 5.1 and 5.2 Excess Pore Pressures at the Thaw Line 
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when, for example, the surface temperature of a homogeneous soil mass 
increases exponentially with time. The curious fact emerges that for 
the constant velocity thaw front, the time factor T decreases with 
increasing time. 

For thawing proportional to the square root of time (n = 1/2), 


the time factor T is constant for all time, and it may be shown to be 


Te, forn= 1/2 (5.8) 


ar® 
where R is the thaw-consolidation ratio defined by equation (3.35) for 
the linear theory described earlier. For n < 1/2, the time factor will 
increase with increasing time. 
For the loading case y' = 0, the settlement curves are very 
similar for S+/Smay less than 0.5. Using the method proposed by Fox 


(1948), it may be shown that these curves are well approximated by 


S 
= == VT (5.9) 
max J 1 

S 

- 50/5 we vee 

max 


For the "self-weight" loading condition (P_ - o4) = 0, no such simple 


fe) 
relationship is valid. 


The excess pore pressures at the thaw line in Figs. 5.1 and 5.2 
are strongly dependent on the value of n. It is clear that in the class 
of problems where n is less than 0.5, excess pore pressure deve lopment 


at earlier times will be more critical, but stability will tend to 
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improve as time progresses. 

The solutions for any combination of the two extreme loading 
conditions may be obtained by simple superposition of the excess pore 
pressures. 

So far, the power law relationship between thaw depth and time 
has been discussed in reasonable detail. Although a relationship of 
this type may be "fitted" to many realistic problems, the function 
itself has no theoretical background in the theory of heat conduction 
in melting soils. A further practical example of the use of the computer 
program in handling an arbitrary movement of the thaw plane is now 
considered. When the surface temperature of a homogeneous semi-infinite 
mass of frozen soil varies sinusoidally, the movement of the thaw plane 
is given by eq. (2.55). It has been shown in Chapter 2 how the 
approximation of the surface temperature sine wave by an average step 
increase in temperature is a reasonable procedure, if the maximum depth 
of thaw at the end of a thaw season is the primary concern. However, 
when one is considering the build-up of pore pressures under such thawing 
conditions, it is the velocity of the thaw front, dX/dt, which controls 
the excess pore pressures to a greater extent. Therefore it is of concern 
to calculate the excess pore pressures in a thawing soil which is subject 
to a sinusoidal application of surface temperature, and compare them 
with the pore pressures maintained in a soil thawing under an average 
step increase in surface temperature. It should also be pointed out 
that in some cases where the thawing of an active layer is taking place 
under natural conditions, it may be more correct to consider the 
surface temperature as varying in a sinusoidal manner. 


Assuming a thaw season of 180 days, and the surface temperature 
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varying sinusoidally between zero and 10°-e; we have 


HEC) 
Se sin( Tag (5.10) 


where t is the time in days. 


Taking the values for thawed conductivity and soil latent heat as 


k 
U 


le 


0.0025 cai/°C.cm.S 


32 cal/em> 
and using eq. (2.55) as derived by Lock et al. (1970), the relationship 


between thaw depth and time is given by 
M162. 4 Sin (=| =OLO7ed Sina! . Sintaae (5.11) 
360 : 360) ° 180 7 


where t is in days, 


and X is in cm. 


This particular example is the same as that used in Section (2.6) 
to determine the depth of thaw under a surface temperature sine wave. 
Equation (5.11), together with its derivative, dX/dt, are used in the 
subroutine "THAW" in the computer program listed in Appendix A.4. The 
functions X and dX/dt, the velocity of the thaw interface, are plotted 
in Fig. 5.5. The computer program prints the excess pore pressures 
and the degree of consolidation settlement as before. The excess pore 
pressures at the thaw line are plotted for two different loading 
conditions in Fig. 5.6, and two different Cy values for each. The 
excess pore pressures rise steadily and reach a peak approximately where 


the slope of the depth of thaw against the square root of time is a 


maximum in Fig. 5.6. 
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k, = 0.0025 cal/*C cms 
lL = 32.0 cal/em3 


Tmax = 10°C 
25 Tg * 0°C 
i: Ste= 0.2261 
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Fig. 5.5 Rate of Thawing for Sinusoidal Surface Temperature 


a 
* 7 i 
i A 
i a rs) 
; ye 
a Acy o 
ee 
re iy Omi 


ee: J Kos 29OR0 © € ie 
, Fina Silos O80) * 4 a 
; Ms at 034 a - on ry r. oy 
* See | ge at Ne 
rase.g a 


af 


te 


| 


ae 
3 Aa8y ‘ : 
SRT RISMaT 4 fw 


ut 7 wags - : 


“DrE.S “f 


u(x,t)/(B + 8X) 


MIALISED EXCESS PORE PRESSURE AT THE THAW LINE 


NOR 


PSVE DU TEM Doi ioe 2 cep uct We Mean) wa yee, 
(R=1) 


0.8 


S=0 


-4 
ae C)= 2.4x10 cm 


sine temp. ace 


O04 


step temp. » 


0.6 
©) 
-4 2 
Cy= 2.4x10 cm/s 
0.4 sine temp. 
step temp. peirateres 
0.2 -4 2 
Cys 9.6x 10 cm/s 
fe) 
Oten2 4 6 8 10 FZ 14 


YM er DAYS” 


Fig. 5.6 Excess Pore Pressures for a Sinsoidal Surface Temperature 
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For comparison purposes, the excess pore pressures obtained by 
thawing under an average step temperature change equivalent to the sinu- 
soidal surface temperature application are also plotted in FUG). 3300. 
The normalised excess pore pressures at the thaw line for a particular 
step temperature application and coefficient of consolidation are 
constant with time, and are shown by the dotted lines in Fig. 5.6. It 
is seen from Fig. 5.5 that the velocity of the thaw front is initially 
greater for the step application of 1s however at later times as the 
Sinusoidal surface temperature rises to its peak, the thaw front 
velocity resulting from the surface sine wave is greater by 
a considerable amount. Consequently, Fig. 5.6 shows that the corresponding 
excess pore pressures at the thaw line rise just slightly above the 
excess pore pressures calculated for the equivalent step increase in 
surface temperature. 

It might be argued that an equivalent step increase in 
temperature will produce a rate of thawing that will always provide 
an upper bound to geotechnical problems such as the build-up and 
maintenance of excess pore water pressures. It could be supposed even 
further that the equivalent step surface temperature would produce 
effects that were much too severe, and would overestimate the adverse 
conditions associated with a fast rate of thawing. It is demonstrated 
here that these arguments are not valid, and that indeed over a 
certain part of the thawing season, the excess pore pressures under 
the sine wave surface temperature exceed slightly those obtained under 
the equivalent step temperature. It is concluded that an average step 
temperature may be used with some confidence when predicting the 


maximum excess pore pressure in a thawing season, in a soil whose 
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surface temperature varies sinusoidally. Nevertheless, if appropriate 
Surface temperature data are available, a more accurate prediction of 
the variation of pore pressure with time is obtained readily using the 


computer program given in Appendix A.4. 


9.3 The Non-Linearity of the Stress-Strain Relationship for the Soil 


Skeleton 

Natural deposits of permafrost exhibit a wide range of void 
ratios in the field. Because the ice matrix in the pore spaces of a 
frozen soil is capable of supporting stress for an indefinite time 
period, a much wider range of void ratios may be expected in the Arctic 
than those encountered in more temperature latitudes. Void ratios in 
fine-grained permafrost soils may range from almost pure ice to a void 
ratio corresponding to the Plastic Limit or less. In many cases, 
extremely ice-rich materials are limited to a fairly small depth below 
the surface. 

However, significant depths of soil may be encountered where the 
soil might still be termed ice-rich. After thawing of these deposits, 
large quantities of water must be expelled from the soil pore spaces 
before appreciable gains in effective stress are to be realised. Stated 
in another way, the constitutive relationship for the skeleton of such 
soils may be quite highly non-linear. 

It is of considerable interest therefore, to examine the 
limitations of applying a linear stress-strain law (such as that used 
in Chapter 3) to soils exhibiting varying degress of non-linearity 
in the soil skeleton. Several methods of incorporating a non-linear 


relation for the soil skeleton might be proposed in an analysis for 
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excess pore pressures and settlements in a thawing soil. Let the 
initial void ratio in the thawed condition before drainage is permitted 
be denoted by Qo? and the final void ratio under the effective over- 
burden loading be denoted by Ce A schematic diagram of a typical 
relationship between void ratio and effective stress is given in 

Fig. 5.7. Also shown are different analytical functions which might 

be used to describe the soil behaviour. Of the four models shown, the 
first two are dealt with immediately, and the remaining models are 
discussed in the next two sections. 

The first relationship between e and o' shown in Fig. 5.7 is 
the linear Terzaghi assumption for the skeleton of a compressible soil, 
and this relationship and its incorporation into the one-dimensional 
theory of thaw-consolidation have been dealt with fully in Chapter 3. 

The second model attempts to simulate the behaviour of the 
soil skeleton to a higher order of accuracy by assuming that the soil 
may expel a quantity of water to obtain a void ratio ey: with no gain 
in effective stress. The skeleton may now consolidate from Qe to 
er until the effective stress is equal to the final overburden stress, 
along a path which is linear but has a slope less than the linear 


model described above. The coefficient of compressibility must now 


be redefined as 


C= 
m =—1_f _._;i_t y (5.12) 
Eee One) eo 


where Mm, is the coefficient used in the linear analysis. 


Rederiving the boundary equation (3.18), for the bi-linear analysis: 
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Fig. 5.7. Different Void Ratio-Effective Stress Relationships 
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] hee 
"Oo TF eym +X ayat (5.13) 


However, the quantity 


e = 
0 oi 


We) m 
1 + eo Mm) 


is a constant, and therefore may be added to the constant eS loading 
term, and the solution may be expressed as a trivial extension to that 
obtained in Chapter 3. For the PS loading condition only (oedometer 


conditions) the solution for the excess pore pressures is 
e -e 
] 
P+ ie marae erf(Rz) 
O tat e,) m 
Tees : (5.14) 
re ae ie eer 
erf(R) + 


vm R 


Although the m, value is now reduced by the ratio given in 
eq. (5.12), and the R value is also reduced, an extra term appears 
in the numerator of eq. (5.14), tending to increase the excess pore 
pressures. This is certainly the simplest method of introducing some 


degree of non-linearity into the solution of the thaw-consolidation 
equation. 
5.4 Influence of a Non-Linear Effective Stress-Strain Relation 


Because high initial void ratios are expected in thawed soils, 


the subsequent stress-strain behaviour of the soil skeleton may be 
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markedly non-linear as shown schematically in Fig. 5.7. Initially, the 
skeleton of such ice-rich, fine-grained soils must expel considerable 
quantities of water with only small corresponding increases in effective 
Stress. As consolidation proceeds under successive increments in stress, 
however, the skeleton becomes progressively less compressible, requiring 
less volumetric strains to gain increased effective stress. As this 
behaviour is thought to be quite common in ice-rich permafrost soils, 

it is of considerable interest to investigate the influence of this 
non-linearity on the consolidation of thawing soils. 

If the soil is thawed and no drainage is permitted, then the 
effective stress in the soil is denoted by ome If drainage is now 
permitted under the influence of subsequent loading increments, then 
the void ratio-effective stress curve may be found experimentally. 
Laboratory tests described in Chapter 4 suggest that for some thawed 
soils tested in this manner, a straight line is obtained by plotting 
€ against the logarithm of effective stress. The same empirical 
relationship was incorporated by Davis and Raymond (1965) when formulating 
a theory of consolidation for normally consolidated unfrozen soils. 

Although the permeability k, and the slope of the stress-strain 
curve Mm both change over a range of loadings, it may be assumed on 
experimental grounds that they change in such a way as to maintain the 
consolidation coefficient Cy approximately constant. Following Davis 
and Raymond (1965), a constant Cy value is adopted together with the 
stress-strain relationship 
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where I. is the slope of the e/log o' line (compression index) 


and o is the current effective stress. 


The general one-dimensional equation of consolidation for a 


saturated soil is then written as 


js bolucn ec id he’oe | nt c Loss 
VileGuesia tt o'} 3x ox Goan ' 


This equation is valid in the thawed soil regardless of loading type or 
drainage conditions at the boundaries. 
For the same loading conditions considered previously, the 


total and pore water stresses respectively are 


CE a ean cyt wy (eel ) 
P =u + ey (See) 


and with the water table at the surface x = 0, the effective stress 


is then the difference between these stress components. 


h = = = ! = 519 
Oo Oo ae iy ry he 0 ( ) 


These equations are valid for any initial effective stress or developed 
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and substitution of these expressions into eq. (5.16) yields 


' ou 
aus, ae ven )a] 
Ut han ellen Pot yx u Ox (5.22) 


A linearizing transformation for this non-linear equation such 
as that used by David and Raymond is only available for the case y' = 0. 
This condition where applied loads dominate will be studied here. 
Numerical means must be employed to solve a loading combination involving 
self-weight stresses. 

Setting y' = 0 in equations (5.19) and (5.22), and applying a 


linearising transformation to the governing eq. (5.22) gives 


awe ow 
ot Vv aye (5.23) 
Pe - U 
where W = 10919 P (5.24) 
The boundary condition at the surface is now 
x = 0; u = 0; w = 0 (5.25) 


At the thaw line, the boundary condition used previously must 
be now rederived for consistency. Previously, a continuity requirement 


generated equation (3.9), that is 


at X°= X(t). (3.9) 


Defining the volumetric compressibility m, as 
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where e is the void ratio, and using the empirical equation 


(5.15) we obtain on differentiation 


0.434 I 
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So, using equations (5.15) and (5.27) it may be shown that 
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(5526) 


(5 2) 


(5.28) 


The boundary condition (3.9) after substitution of (5.19) and 


(5.28) becomes 
atexXt = X(t); 


ou 
1 + RT aU Ce oye SU ’ Cy ax wee 
0458 has 2210 a’ 
°6 dt 


(5.29) 


Setting y' = 0, provides the boundary condition required to 


study the applied loading condition. Rewriting in terms of w(x,t), 


the transformed variable: 
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(5.23), (5.25), and (5.30) in terms of the transformed variable w 


are entirely analogous to the corresponding equations derived for the 
linear theory in terms of u(x,t). Consequently the solution for the 
linear thaw-consolidation equation may be applied directly to obtain 


w(x,t). For /t thawing 


y= wv (z,t) . erf (R z) 


‘sale Z 
Oi. -R (5.31) 
10919 (F2 } erf(R) + © 
7 R 
where z= Wty 
(5.82) 


and the excess pore pressure is obtained from the modification 


ca ARE (%) 


ide wor, (5.33) 


Equations (5.31) and (5.33) form the exact solution for a 
thawing soil with consolidation occurring solely due to an applied 
load Po: The soil has an initial void ratio in the thawed state 
corresponding to an initial effect stress of Oy * 

The effects of self-weight stresses may be determined relatively 
simply by numerical methods, and it is considered sufficient for the 
scope of this study to examine the effects of an applied load upon 
the excess pore pressures. 

The ratio P cho is a measure of the effective stress change 
in the thawing soil, and therefore of the non-linearity of the effective 
stress-strain curve. This ratio now enters as an additional dimensionless 


quantity, and its effect on the excess pore pressure distribution in the 
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thawed soil may be investigated. As in the linear theory, the normalised 
excess pore presure distribution is unchanged throughout the thaw- 
consolidation process, although it is now dependent on the stress ratio 
P/F: The degree of settlement in the thawed soil however is 
independent of time and the ratio Eo This fact was also recognized 
by Davis and Raymond (1965), for this type of stress-strain relationship. 
Consequently, the linear theory may be used equally well for calculating 
the degree of consolidation settlement. 

Figure (5.8) shows the excess pore pressure at the thaw front, 
as a function of the stress ratio and the thaw-consolidation ratio R. 
It is seen that as the stress increment ratio increases, the excess pore 
pressure condition at the thaw line increases to a more critical level. 

In cases where high initial void ratios are present in the 
soil on thawing, the excess pore pressures may be underestimated 
considerably, if only a linear relation between void ratio and effective 
stress is considered. The logarithmic function given by eq. (5.15) 
provides a method of accounting for any required degree of non-linearity. 
The function however always requires some value of Opt the residual 
stress. In very ice rich soils where this extension is obviously most 
useful, ome is immeasureably small or zero, and so apparently cannot be 
used in the analysis, as the stress ratio is now infinite. However, it 
must be remembered that the reason for the introduction of such a 
function shown by the dotted line in Fig. 5.7 is to model in some way 
the high degree of curvature between ey and some void ratio e,- If a 
sample of frozen soil is thawed under a small nominal loading (often 
the weight of the load cap in an oedometer), it consolidates to a void 


ratio e, shown in Fig. 5.9. Now the stress path taken by the soil to 


at 


thie ts 


shila ars aca wean ault at ah : 
~werld odd suoruars bopnertomy - ae ud Pita’ > 
ates 2esit2 ont no Siabnsyeb won at a dp 0 jis 
ai saver ftoz bowsrt ond ar irom | e, “ 
bastnpoosy o#fs 2ew Jost eNiT .' ol, 4 fore at bas aes to snsbna 
_giienottetey mbete-eeevte to a0Ns erat oF (¢aaen), bre 
onitetuotso, yot flow MET supe beau ad yon yNoss 48 ' 
. Sneniat tase not Sab Tozne> Yo notes 
,inot? wed att ts swezarg syeq 2¢a9xo ond wore (8. é) ugh 
AW ottey notsebtToenos-wans ant bre. Oft6a caente, ond Yo: sottomt 
ay0q 2e99xs ed} , eozsaiant often SH omg NT ezonse ‘rts 26 des soe 2 a 


a 


/fevs : featding ayom 6 Os eaessroni oni wai ad3- $a Hots thn « 


panauraeeienan od (ysm coweesg or: cn D9 
avtsostis brs often brov no4wited notes lay ; 7 iit 

(at.2) pa yd hevie norton ott 

\“ rrsontl-nat to Sevpob beriuper a] mot el 


faubteay oft 4" 0 Fo oulay ome, =a 
szom yt 2 uel ¥do 2t notenatKe ef - sity ef de 
ad Jorg) Vansasggs 08! brs ones ad ailcncsenei ee fu 
+} ewevewotl “cotta? won at otter ative ott ee, 2tey i sits bit 
borgat Bat wot ages anit: $eeit bervodemsry 9 
ys omoz nt f abort of at % ‘ sshieidecinihaae ant) ed manele 


Pi 


571. “et otter bioy sno2 big lk nema ows srs to sone fot a ond a 
masto) een ik sine’ bos ww bewed: at ee sto 


B fave: to m6. “f5) 


1?3 


AAOBYL ABPBULL-UON UOJ SUL] MEY] 7e SauNSSaug B0g Q's “BLY 


Y OILVY NOILVGITOSNOD -MVHI 
Ol G 0'| SO lO 


N 
oe) 


A 
o) 
SYNSS3Yd 3YOd SS39X3F 


(O-ayn 


= 
fo) 


80 


JNIT MVHL 3HL Lv 


01 


E1@* ag OM pleeciee Sp dUSM FyUE Lok yOU-jiuEsL jyEOKA 


BViIjO 


LHVAM- COMZOTIDVLION 


iQ: 


ro 


Oe 


174 


arrive at this void ratio is unknown, but it does not appear important 
to know it in any case. Additional load increments are now placed on 
the sample, and the void ratio at the end of consolidation of each 
increment is determined in the usual manner. A typical set of data is 
shown plotted both on natural scales and logarithmic scales in Fig. 5.9. 
If the points plotted on the e vs. log o' diagram are now extrapolated 
back to meet a horizontal line drawn through Qo» the initial void ratio, 
a complete functional relationship between e and o' is obtained. This 
construction is shown by a dotted line in the e vs. log o' plot in 
Fig. 5.9. A very small value of Oe is now obtained, and whether or 
not this value is the actual value for the residual stress in the 
undrained thawed soil is probably immaterial, as we are solely 
interested in obtaining a value of the stress ratio PhO, that will 
express the high degree of non-linearity associated with the soil 
compressibility behaviour. The "fictional" portion of the e vs. 
log o' curve between ey and e, is also plotted on the graph of e vs. 
o' ona natural scale. It is seen that for practical situations the 
fictional value of on is so small that the first part of the e vs. 
o' curve is approximated to a high degree of accuracy. The excess 
pore pressures in a thawing soil exhibiting this e vs. o' relationship 
would be determined by using a stress ratio equal to the final effective 
stress ae divided by the value of Ty's which in the case shown here 
is 0.02 kg/cm. 

It is concluded then that the use of a logarithmic variation 
of effective stress with void ratio may be used with considerable 
advantage in the theory of consolidation for thawing soils. When the 


soil is initially at a high void ratio, and subsequent consolidation 
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causes large strains in the soil skeleton, a residual stress that 
may have little physical basis may be determined. This enables the 
soil stress-strain behaviour in the initial highly compressible 


condition to be modelled mathematically. 


9.5 The Thawing of Ice-Rich Incompressible Soils 


The fourth model in Fig. 5.7 relating void ratio to effective 
Stress is visualised as a bi-linear relationship which is vertical 
initially, and on reaching a final void ratio Crs proceeds horizontally 
to the final effective stress. In many cases of practical interest, 
the solid constituents of the frozen soil are of coarser grained 
material. Thawing of such soils permits the liberation of any moisture 
which is in excess of the final void ratio er. Further loading of a 
coarser-grained soil produces little additional compression, and so 
consideration of a model such as this provides a useful insight into 
the thawing behaviour of ice-rich incompressible soils, and also gives 
a convenient lower bound to the series of void ratio-effective stress 
relationships which it is possible to consider for any soil. One 
might envisage a further application of this stress-strain relation 
involving a finer-grained soil. 

In a clayey silt, for example, it is common to find rythmic 
ice banding in the frozen material. If the interbedded soil layers 
are highly overconsolidated with respect to the present overburden 
stress, and would themselves demonstrate little compression on thawing, 
then perhaps an approximate analysis could be carried out by assuming 
a rigid incompressible soil skeleton accompanied by uniformly 


distributed excess ice. However, it is clear that upon adoption of 
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this model, the settlement ratio S,/S is always 100%, provided 


max 
that the discharge capacity of the soil is not exceeded (i.e. the 
effective stress at the thaw line does not become zero.) 

The excess ice content of an ice-rich incompressible soil is 
defined on the basis of thawed void ratios as follows, 


(2 Sk 


Ere oe 
6) 


As the soil skeleton in the thawed zone is considered 
incompressible, the equation governing dissipation of the excess pore 


pressures in that zone is the Laplace equation in one-dimension. 


du 
<5 = 0 (5-35) 
dx 
The solution to (5.35) is written simply in the form 
u = Ax + B (5-36) 


At the surface, a free drainage boundary exists, and so 


(5°37) 


w 
i] 
Oo 


At the thaw line, the continuity equation (3.9) derived previously 


is written down 


k Ou 
AV aX 

Ato =ex(t)s Sia, Gaal (3.9) 
Mer 


As the soil is incompressible and therefore possesses no capability 
for transient storage of excess pore fluids, all excess water must 


be expelled from the thaw front, and so 
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a eae coed (5.38) 


and substituting eq. (5.38) into (3.9) the boundary condition at the 
thaw line is obtained: 
DU ee ttt cae ex 


at x = X(t); are " at (5.39) 


Combining eq. (5.36), (5.37) and (5.39) the solution is 


u= 41 a . (5.40) 
Normalising the excess pore pressure and the depth x, the solution is 
rewritten as 


Yuy EY g dx 


pee Ue = —— 
Po+yXx “k(P ty x) * dt Sy, 
fe) fe) 
It is desirable to compare the results obtained from this incompressible 
analysis with those obtained using the linear theory described in 
Chapter 3. This is only possible if we define an equivalent coefficient 


of volume compressibility for the incompressible analysis to be 
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The inverse of this term appears in eq. (5.41), and therefore a method 
of comparison with the linear theory of thaw-consolidation suggests 
itself. 


For thawing proportional to the square root of time, 


X dX = oe 
dt 2 (5.44) 

and at the thaw line 
Zea | (5.45) 


and so, placing the equations (5.43), (5.44) and (5.45) in eq. (5.41) 


provides 


(5.46) 


To facilitate the comparison, the analagous thaw-consolidation ratio 


is defined as 


seh 


Va 


and from eq. (5.46) the much simplified expression is: 


* 
ey = 2 2 (5.47) 
0 iY 
Equation (5.47) allows the expression of the normalised excess 
pore pressures at the thaw line in terms of an equivalent R value, 


which applies to an incompressible soil with excess ice. The relationship 
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is shown plotted in Fig. (5.10) together with the corresponding 
relationship for the linear theory of thaw-consolidation. The 
relationship expressed by eq. (5.47) is independent of the self-weight 
ratio We At the lower R values, the excess pore pressures predicted 
by the two theories are very similar. However, in the middle range 
of R values, the incompressible model predicts a marked increase, and 
at R = 0.7 the excess pore pressures become equal to the effective 
overburden load. It must be remarked that although the incompressible 
model appears to predict pore pressures that are independent of the 
self-weight ratio, the m value and hence the ae value contain the 
quanti ty ee + y'X). This term varies with time if y' is not zero, 
and so therefore does the R- value. 

Under conditions of no external loading ee = 0), the 
discharge capacity is exceeded in all cases to a certain finite depth. 
As the discharge capacity is proportional to the quantity y'X, at 
small depths of thaw the rate of liberation of excess water exceeds 
the small discharge capacity. The depth to which the thawing soil will 
remain totally unconsolidated may be found by establishing an inequality 


from eq. (5.47). That is 


u *y2 5.48 
pty = a(R )S > (5.48) 
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or 

D * 

m 
coma aa TaLe VN (5.49) 
ete aaK 


* ° ° 
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2 
oe E. 
ee (5.50) 
or 
af VReLe 
Yee Oy 


ék y" (5.51) 


If the depth of thaw is less than this critical depth of X, 
then the thawing soil is totally underconsolidated for the self-weight 
(y') loading conditions. This is the most striking difference between 
the incompressible model and the results predicted by the linearly 
compressible relationship, which always predicts a constant degree 
of consolidation at the thaw line, with no critical depth of thaw such 
as that given by eq. (5.51). This preliminary study would suggest that 
for soils exhibiting a high degree of curvature of the stress-strain 
curve, and which are thawing under self-weight conditions, excess pore 
pressure conditions would be most critical at early times, and would 
improve with depth. 

The excess pore pressure conditions can be greatly improved 
in the analysis described above by the simple provision of a surcharge, 
Po: By the same argument developed above, the placement of a surcharge 
of magnitude 

ae E. 


Ps = 5 alread) (5.52) 


will increase the discharge capacity to the extent where some degree of 
consolidation is obtained from the initiation of thawing onwards. 
In summary, the stability of thawing ice-rich soils which have 


an incompressible soil matrix may be considerably improved by the 
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provision of a surcharge loading. This has the effect of increasing 

the discharge capability of the soil to an acceptable level. If, 
however, the surcharge loading were not distributed uniformly on the 
surface, as in the case of a berm, then additional shear stresses would 
be induced in the subgrade. Such shear stresses might reduce the 
stability of the subgrade, and thus negate some of the beneficial effects 


of the surcharge loading. 


9.6 The Consolidation of a Thawing Soil With a Layered Profile 


So far the solutions to problems in the one-dimensional thawing 
of permafrost soils have confined themselves to the treatment of a 
Single compressible layer exibiting a homogeneous set of geotechnical 
and thermal properties. In many cases of practical importance, the 
depositional and thermal history associated with a natural deposit 
of permafrost may be such that extreme changes in consolidation and 
thermal properties occur with increasing depth. These changes in 
properties may occur gradually, but are more likely to be quite 
discrete. A sudden change in depositional history for example may 
bring about a distinct change in soil type. In such instances, two or 
more discrete layers of soil may be assumed present, each possessing a 
different water content, latent heat, permeability and coefficient 
of consolidation. 

It has already been shown in section 2.7 how the movement of 
the thaw plane through a soil profile displaying two different sets of 
thermal properties is analysed approximately. Once the history of 
the thaw interface is established, attention may be focused on the 


associated problem in soil consolidation. Figure 5.11 shows 
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Fig. 5.11 Thaw-consolidation in a Two-layer Profile 
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schematically the consolidation of a two-layer profile, each layer 
having the properties Cys K, k and L, the coefficient of consolidation, 
permeability, thermal conductivity and the volumetric latent heat 
respectively. 

The movement of the thaw interface X(t) through the two-layer 
profile is derived in terms of the Stefan formulation in section 2.7 


and is rewritten here as 


Mee I, 
X =o/t = — Vt 3) X<h (2.60) 
1 


where ue is the suddenly applied increase in surface temperature. 


On penetrating the second layer, the motion of the thaw plane is 


described by 


X >H (2.63) 


where ty is the time to thaw the first layer of thickness H, 


and is given by 


BL ic eee 
ty = H/a 


The consolidation problem is summarised formally as follows for t > toe 


when the thaw interface has entered the lower layer. 


In layer 1, 
mt ay 0 H (5053) 
t > ave NEL ve me << : 


t 220: Daas Oi x = 0 (5.54) 
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In layer 2, 
; , dU5 au, 
ates =z UFC : Hest k(t) (5.55) 
fe) ,) Vo ayo 
ta 
V5 OX 
ca t)3 gp - oo + eX! = U, x ee KC) (5.56) 
dt 


At the interface between the two layers, the excess pore 


pressures are equal and so 


Uy = Uy3 x =H (5.57) 


The equation of continuity for excess pore fluids must also be 


satisfied giving 


te iy K x =H (5.58) 


0 10x 


The system of equations (5.53) to (5.58) describe the dissipation of 
excess pore pressures in the two layer soil mass for X >H. The 
solution to the consolidation problem for X < H is known analytically 
from the results presented in Chapter 3. The equations presented above 
can be generalised to treat the problem of thawing in a multi-layered 
system. If the thaw plane is passing through layer n of a multi-layered 
orofile, there will be (n-1) equations of the type expressed by 
equations (5.57) and (5.58) at the (n-1) internal interfaces. The 

same equations for the free drainage surface at x = 0, and for the thaw 
plane at x = X(t) will apply as before. It is only necessary, however, 


for the scope of the present study to examine some of the effects of 
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the passage of a thaw plane through a two-layer system. 

The equations (5.53) and (5.58) are written in finite 
difference form in preparation for solution. The details of the 
numerical method are given in Appendix A.5. A computer program was 
written to solve the resulting finite difference equations for the 
excess pore pressures u, (x, t) and un(x,t). The settlements are computed 
as usual by numerical integration of the excess pore pressures using 
the trapezoidal rule. A listing of the computer program and the method 
of data entry are given in Appendix A.6. 

The program is used to investigate the effects of a surface 
layer having different consolidation properties on the maintenace of 
excess pore pressures in the underlying thawing layer. Two simple 
cases are solved for an applied surface loading where the upper layer 
has a lower coefficient of consolidation and permeability, and also 
where the coefficients of the upper layer are considerably higher. The 
same thermal properties are assumed to exist in each layer, in order to 
examine the effect of a change in geotechnical properties alone. The 
excess pore pressures while the thaw plane is in the upper layer are 
calculated from the linear theory in Chapter 3. The excess pore 
pressure profiles are shown in Fig. 5.12 and Fig. 5.13. for typical 
properties exhibited by "silt over clay", and "clay over silt" 
situations. As the thaw plane enters the underlying layer, a sudden 
change in the excess pore pressure at the thaw line occurs. Each 
figure also denotes by an arrow the excess pore pressure which would 
be obtained if only the lower layer were present. 


The excess pore pressure in the two-layer profiles apparently 
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approach this value in an asymptotic fashion, but in fact after some 
years of thawing, the effects of the surface layer on the thaw line 
pore pressures are still pronounced. The conclusion migit be drawn 
from the results of these two simple probienetenee the prescence of 
a surface layer which participates in thawing and the subsequent 
consolidation, although possibly minor in extent, may have long term 
effects of the thawing of the underlying layer. 

A more realistic field situation is considered in order to 
demonstrate the applicability of the computer procedure in Appendix A.6 
to the solution of practical problems. The excess pore pressures are 
obtained for the thawing of an ice-rich silty clay soil which is 
overlain by a five foot layer of peat. A surface temperature of +10°C 
is applied to the surface of the peat, which has widely differing 
thermal and consolidation properties from the underlying clay stratum. 
The rate of thaw and the thaw line pore pressures are computed also 
for the case where the peat layer is removed, and the same surface 
temperature is applied. In this simple manner, some of the effects 
of the removal of a surface organic layer on the stability of a 
thawing soil may be demonstrated theoretically. The thermal and 


consolidation properties for the two layers are assumed to be as 


follows. 
Peat 
Water content w = 550% 
Submerged density —y'_ = 3.5 1b/ft® 
Depth H SO feet 
Latent heat L, ell cal/em> 
Thawed conductivity k, = 9 x 107*ca1/°C.cm.S. 
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Surface temperature Ue = bade 
Thaw rate in peat alone Quetac oe ft/year |/@ 
Consolidation coefficient Coe 2200 ft?/year 
Permeability K = 130 ft/year 
Silty clay 
Water content w = 50% 
Submerged density y' = 45 1b/ft® 
Latent heat L, = 46 cal/em® 
Thawed conductivity k= 27 x 107" cal/°C.cm.s. 


Thaw rate in clay 
alone (if peat removed) a = 6.32 ft/yr'/@ 


Consolidation coefficient Cc = 36 ft2/yr 


iT] 


Permeability K 1+0-fFt/yr- 


The excess pore pressure profiles at different times for 
thawing in the two layer profile are shown in Fig. 5.14. Also shown are 
the thaw depth with time and the normalised excess pore pressure at 
the thaw line. The same relationships are shown for the case where 
the peat layer is removed, given by a dotted line in Fig. 5.14. After 
ten years of thawing, it is seen that the thaw penetration into the 
silty clay layer has been reduced by some fifty percent. More 
significantly, the build-up of excess pore pressures at the thaw line is 
such that only about fifty percent of the pressure is attained of 
that value which would be calculated if no peat layer were present. 
The prescence of the peat layer is, in general, responsible for three 


effects, all of which serve to make less critical the stability of 


ae ba ssi i at 


nade " hous» re 
sot eer = ; 


YOe ey 
peat By 
“St\di @h = 'y 

Saati es. ah = gt 
12g Dei ? BE xO ¥S 


% 
Ma 


SMa? 96.8 
ahah at 


Sia 0. ; 


a 
2. .6 


tt 
= 


ais wore oatA ar, a ae “nt ee | 
ts SsiN2eagng) arog. ak 
svattin 9289 okt Yo? mwore . sha Te 
ASttA- BE.d iy i, ant bad ab | Cao 
atis oink not tevsens ale a a by 
ato .dHias18g eo Kiana cmt eat ns e 

2t sitl wart. ong +6: + eae 2aepe. me 


bors 10% serssuii 5 dela - yous | 
Bo vit lidst2 eid, Iaptoine 2eat gam 


193] 


U 
NORMALISED EXCESS PORE PRESSURE BH+8X 
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Fig. 5.14 Excess Pore Pressures for 'Peat Over Clay' 
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the underlying thawing layer, 

(i) A high latent heat barrier is available to impede 
thawing of the upper layer. 

(ii) After thawing of the peat, its lower thermal conductivity 
reduces the rate of thaw in the underlying stratum, so reducing the 
overall thaw penetration, and the excess pore pressure maintained in 
the underlying layer. 

(iii) A relatively free draining boundary is provided, and 
negligible excess pore pressures are maintained at the base of the 
peat layer. 

If the peat were removed and replaced for example by a layer 
of gravel, conditions (ii) and (iii) would still be available to 
improve foundation conditions, however the latent heat barrier 


described by (i) would in all probability not be available. 


9./ The Analysis of a Compressible Soil with Discrete Ice Lenses 


When a fine-grained soil freezes, a combination of the 
ambient conditions may prove extremely favourable for the formation 
of segregated ice. Some of the conditions favourable for the creation 
of ice features in frozen soil have been discussed in Chapter 1, and 
they include the slow rate of heat removal from the freezing plane, 
the availability of water to the freezing plane, and a low overburden 
pressure at the level of freezing. The conditions for the formation 
of segregated ice near the surface of a deposit of a freezing soil 
are very often realised. In soils with silt-sized particles, 
segregated ice in the form of discrete layers or bands often occurs. 


In the finer grained plastic soils, ice banding is also found, with 
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the bands oriented perpendicularly to the maximum principal stress. 
However, vertical ice veins and other ice forms are also found regularly 
in the soils with more clay sized particles. 

So far a complete analysis of the position and magnitude of 
ice layers or even the amount of distributed ice which may be found in 
a freezing soil has not been carried out. It would therefore seem that 
the analogous problem in thawing soils would defy analysis as well. 

The interaction between the heat and mass balances at the frost 
interface, the transfer of heat in the frozen and thawed zones combined 
with the thermodynamic relationships for ice/water/soil interfaces seem 
to preclude the use of a simple mechanical description of the problem 
of freezing in fine-grained soils. Hardly content with the problems 
involved with these considerations, some authors proceed even further 
to attempt to describe the moisture transfer in unsaturated soils. 

The solution of problems involving the thawing of soils 
containing horizontal ice bands is rendered considerably simpler by 
knowledge of the position and extent of the segregated ice forms 
initially. A situation involving a layer of saturated compressible 
soil overlying a single layer of ice will be analysed. It is obviously 
impractical to suggest that the analysis developed here be used 
generally to predict the excess pore pressures and thus the stability 
in a thawing soil containing ice lenses, but hopefully all the 
important variables might be contained in a dimensionless statement 
which would indicate the potential stability above a thawing ice layer. 
Such a dimensionless statement would not only provide a stability 
index such as the 'R' value derived in Chapter 3, but should also 


suggest remedial measures which might be applied to improve Stability. 
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At the outset, a solution is considered for excess pore pressures 
in a thawing soil over a layer of ice which extends downwards indefinitely. 
This assumption is made for simplicity at this point, and it will be 
seen later how the procedure is modified to include any ice layer of 
finite thickness. It is assumed that a step increase in temperature has 
been applied at the surface and that the thaw plane is moving through 


the frozen soil at a rate given by the simple Stefan formula 


yt (5459) 


whe re Le is the volumetric latent heat of the soil. 


The excess pore pressures obtained while the thaw plane is still in the 
soil are given by the linear theory of thaw-consolidation in Chapter 3. 
For the case of self-weight loading only, the excess pore pressures 


exhibit a linear profile with depth, and can therefore be written as 
u = B y'x (5760) 


where B = —;— and is a constant with time. 


When the thaw plane encounters the ice layer, the height of 
soil solids above the thaw line remains constant with time. If we 
now assume that the strains in the soil skeleton associated with any 
increase in excess pore water pressure (caused by the extra influx 
of water into the soil) are small, then the distance between the 
soil surface and the thaw plane becomes constant. This assumption is 


expressed schematically in Fig. 5.15, and essentially implies that 


samizeayg eAog ple bee bovebiarigg: pF ‘ae! acai 
oufas FAT PRAT er | ‘gonisdxs dotew sof Fo" ae ‘ 4 4 
od Fit St Bas. dntog 2tdt 36 wtottgate % r o mal 

io vayet aot Yas sbulane of battibon’ cae sive . eo 


ze oyudsysamed nt ‘senonont= gate’ 8 fade boride 41 ean 
dguortd onivon ai bake world afd. ede Ben ‘goptwe - 


26 hotdiw 9d eee go) Bete ‘ate, 2 amcatTe wsonit bl ate 
(08.2) pirate rey c- 


ois gti dies ates  abibie ah a 9180 


anal 
aia ‘ails ag wad tt nan ‘Carey 
ni sie vSrehT wes, wrtd evods eb i Tae Tod 
vite ee Pr bets hoses. noget ave. Tee afin antents anid teeta smuees wort 
Mtn? adye Sit ya boeuso) sweeetg wstaw Bog seams at seasiont 
eas neoniod sonsietb: oft. oath fone ore {Phoe wes oder nadew ae 

ai. sotiqmiees 2i¢ sSMEIENGS eased, ret wert sed bs sostwue Phoe 
~ Beds died vf fersng229) brs’ Be ii ee, ‘beaeetaxe 


totep ted “sry oN aah ipsa ssa 


ow TE. anh) Adiw feREmD on 2 


195 


Oxia 
> og hit 


Oc Ne a) 


NENENS 


pb § 
mas 


» 


= 


=\ 


b) Thaw plane in ice 


Thaw plane in soil 


q:) 


Fig. Sc 15 “thawing in a Soil-Ice Profile 
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the soil is not absorbing an excessive percentage of the water 
discharge from the thaw plane to the surface. It will be seen later 
from the solution to the problem that this is a reasonable assumption. 
If the distance Xo from the surface to the thaw plane is now a constant, 
then certainly a steady state temperature distribution exists in the 
taawed zone, where the temperature is Ry at the soil surface and Te 

at the thaw line. Applying the usual heat balance equation at the 


thaw line: 


[a a (5.61) 


where an is the rate of advance of the thaw plane into 
the ice layer, 


L is the volumetric latent heat of water, 


W 
ane Aig ett.) (5.62) 
DON ests =F 
Vv 
OX Ao 
From this, the rate of degradation of the ice layer is 
dX = Eaels > Te) 
t Li a (5.63) 


This is a constant velocity with time, as distinct from the solution 
to the heat transfer problem through soil, where the velocity dX/dt 

is proportional to time !/é, It is clear that the thaw front velocity 
will be reduced when an ice layer is encountered, as be is greater 
than Los the latent heat of the frozen soil. At first sight this 
might be thought to improve the excess pore pressure conditions in 


the soil above the ice lens, however it must be remembered that all the 
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melt water from an ice lens enters the soil as excess pore fluid, 
whereas when thawing through a mixture of soil and ice, only a portion 
of the melt water may be thought of as excess pore fluid. 

As no new thawed soil is entering the thawed zone, the usual 
boundary condition expressing the continuity of water at the thaw 
line is no longer valid, and the problem in soil consolidation concerns 
a fixed height of soil Xo: At the thaw boundary at the surface of 
the ice lens, the continuity of the mass of water is satisfied as 


follows: 


Y Se KE: t>t (5.64) 


where to = x70, and is the time to thaw the soil above the ice, 
K is the permeability of the thawed soil, 


and Y; is the unit weight of ice. 


Substituting the rate of melting of the ice given by eq. (5.63), we 


obtain 
UK Kt — = (5.65) 


Letting z = = (5.66) 


0 
and dividing eq. (5.65) by pe + 7 to normalise the excess pore 


pressures 


) u/(P. a y'Xq) Vou iTg - Te) 


0 Z K(P + a) L (5.67) 
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The right hand side of eq. (5.67) is a dimensionless quantity 
which shall be denoted by 


nee Yj oe # T.) 


“RF YEY T (5.68) 


W 
The dimensionless excess pore pressure gradient applied at 
the base of the soil after encountering the ice lens is given by 


eet ol aA (5.69) 


The solution for the excess pore pressure u(z,t) above a 
thawing ice layer may now be found. The linear Terzaghi consolidation 


equation is assumed to govern the dissipation of excess pore pressures, 


therefore: 
r) 32 
QO: 247404 hetero 4k Tcxd0 (5.70) 
3 2 
OZ 
cyt 
where: oT =) aa 
a4 
0 


and the time ty is set equal to zero for simplicity. 


Within the assumptions of the Terzaghi theory, the coefficient of 
consolidation is an independent constant, and so the value of Cy, 
used here is the same as that value used when the thaw plane is 


in the soil. The surface is free-draining, and 
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At the base of the fixed height of soil, a constant influx of water 


provides the boundary condition (5.69) 


ou 


= Y 1 = MK 
7 adel Diy ka ae A; T> 0 (Sa 72) 


The case of self-weight loading will be considered from now onwards, 
due to the simplicity of the initial values. The applied loading 


condition must necessarily be analysed numerically. From eq. (5.60) 


the initial values are written as 


UF 
(ee) 


OR az 5 = BZ J 


(5.73) 
Y Xo 


As the equation of consolidation is linear, and the initial values 

form a linear distribution, the solution is simplified by taking 

the initial values from eq. (5.73) as a new reference line for the 

pore pressures, and then consider the further disturbance to the pore 
pressure distribution by the normalised gradient A applied at the base. 


This is performed mathematically by making the substitution 


= U - 
W y'X a4 (5574) 
1@) 
and 0) ea Le 
Oz aon 9z B oe 


The surface condition (5.71) and the governing equation (5.70) remain 


unchanged in terms of the transformed dependent variable wet). The 


boundary condition at Xo becomes 
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2 =: One (A - B); Te .0 (5:76) 


The initial values become 


Oeenci cris es Oe Nee co (5.77) 


The solution to eq. (5.70) subject to equations (5.71), (5.76) and 


(5.77) is obtained from Carslaw and Jaeger, 1947, page 113, and is 


2 
w(z,t) = (A-B)dz- D> 2 sin ze" | (5.78) 


where M = (2n + 1)1/2 


Replacing w by u(z,t), the excess pore pressures are 


foo) 


2 
Eee ae wenn (A = a aD bs a Sin(M ze" e (5.79) 
M 


ox 
0 n= 0 


=e 


where A andB are the dimensionless physical constants 

defined by equations (5.60) and (5.68), and which 
express the initial and final excess pore pressure gradients in the 
thawed soil. This solution converges slowly for small values of the 
time factor T. The solution (5.79) was evaluated by summing the terms 
in the infinite series on the computer. Figure 5.16 shows the 
normalised excess pore pressure at the thaw line plotted with the 


time factor for different A values. The value of the initial pore 


(ata) GT 


bas (ave). w(t. a} 


(eta) 


edass anos fasted 28. | 

dotdw. brs (88, 2) ‘brs (08.8), amottu 
ad At 2tnsibsng. sweestg 10g, 2289) 
eid t6 zeut sy rhame vor a 28 a 
ai St ond. shit vd ce a ve 


fe) ce 
= x S 


oxas(EX}N BSAVT SOI LV 


1.0 


10" 
2 


TIME FACTOR T= Ct/X, 


T ss) ro) 


oe) ro) 
3S3YNSS3SY¥d 380d SS3IDXK3 


201 


Fig. 5.16 Pore Pressures at a Soil-Ice Interface 
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pressure gradient B is taken to be zero, as it is desirable to isolate 
the effects of the thawing of the ice layer alone. The corresponding 
curves for any other positive value of B are easily evaluated from 

eq. (5.79), or they may also be obtained from curves for B = 0 given 


in Fig. 5.16 by the use of the relation 


u = A- ("52) A-u 
B*0 


B= 0 (5.80) 


where su represents the normalised pore pressure for 
= 0 
B = 0 given in Fig. 5.16, 
and eu represents the required pore pressure for any 


BO 


non-zero value of B. 


The preceding solution is valid where the ice layer continues 
indefinitely below the soil. 

When the ice layer is of finite thickness, as is usually the 
case, the time to thaw the layer completely is calculated from eq. 


(5,63) as 


Poe eee (5.81) 


where h. is the original thickness of the ice layer, 


and t, is the time required to thaw the ice layer completely. 


fi 


This time may be converted into a time factor by 
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Knowing the time factor at the conclusion of thawing of the ice lens, 
Fig. 5.16 may be used to determine the final (and usually the worst) 
pore pressure conditions over the ice layer. 

It is seen from Fig. 5.16 that if the value of A is less than 
unity, then the normalised excess pore pressure will rise to this 
value, and maintain in the limit a pore pressure numerically equal to 
A. If, on the other hand, the A value is greater than unity, then 
stability is only assured for a finite period of time, and eventually 
the excess pore pressures become equal to the effective overburden 
weight, and complete instability results. 

It appears then that stability of a compressible soil mass may 
be gained in two ways. If the A value is less than unity, then some 
effective stress is assured indefinitely at the thaw line. When A is 
greater than unity, instability results unless by calculating Te 
from eq. (5.82) it may be shown using Fig. 5.16 that the normalised 
pore pressures never attain their maximum value. It is seen that Te 
is proportional directly to hss and inversely to Xo: This indicates 
that only small ice layers may be tolerated close to the surface, and 
that larger layers might be thawed quite safely at larger depths. 

The presence of the term De u Ye Key on the denominator also 
reinforces this statement concerning the increase of stability with 
depth. This also suggests a method of increasing stability. If a 
surcharge loading ae were placed over the thawing ground containing 
ice lenses, the A value is reduced. If large ice layers are present 
near the surface, and instability is predicted at low values of Coke g 


the placement of gravel fill might well be used to reduce the A value 
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below unity, and thus render a foundation considerably more stable 
on thawing. 

In order to enhance the explanation of the relative effects 
of various parameters on the stability of the ice layer, a samp le 
problem is solved in real dimensions. To lend some 
realism to the sample solution, a set of soil properties shall be 
adopted which are the same as those found at the Inuvik test pipeline 
site (see Chapter 6). This soil is a clayey silt, and the following 


set of representative soil properties are used 
5 


permeability Kv eee 2 eoexe0e wens 
coefficient of consolidation c= 1.1 x 10° om?/s 
water content w = 40% 

Submerged unit weight y Ta Usoe g/cm? 

soil latent heat L. = 41.4 cal/em? 

thawed conductivity Rene ine cal/°C.cm.s 


Two rates of thawing shall be examined; one corresponding 
to the fast rate of thawing under a hot pipeline, say; 

surface temperature Te Sava hed od 
and the other corresponding to a natural rate of thaw which might take 
place in the active layer; 

surface temperature dies = 12%. 
Using the Stefan solution, the a values for thawing in the soil are 
calculated to be 


a, = 0,093 cm/s!/2 for the high T, 


and Qo = 0.038 cous for the low lite 


Combining these with the Cy value provides the two R values: 
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Ry = 0.443 


¥ 


Boe Su Ost 
The pore pressure maintained in the soil, and therefore the initial 
values for the ice layer problem are 

B ee 0 283 


B 0.061 


2 
The effects of an ice layer in this thawing soil are now considered. 


For each rate of thawing cited above, the top of the ice layer is 
assumed to be at Im. or at 3m. For Xo = Im, the A values for the 

two rates of thawing are 1.0 and 0.169 respectively. The relationship 
between time factor and time may be calculated from eq. (5.70) for 
each depth of ice layer, nee Figure 5.16 and eq. (5.80) are now 

used to plot the excess pore pressures at the ice/soil interface with 
time. Figure 5.17 shows the normalised pore pressure plotted with 
time in days for the two thaw rates and depth of ice layers. The 
history of the thaw interface is also plotted for each case. 

For the fast rate of thawing, the pore pressure above an ice 
layer at a depth of Im. rises to 94% of the maximum effective stress 
within 10 days. This represents a thickness of 20 cm. of ice thawed 
in this time. An ice layer at 3m. depth under these thawing conditions 
produces pore pressures that are not significantly more critical 
than those produced by thawing in the soil. 

For the slow rate of thawing, the prescence of an ice layer 
at Im. depth does not produce pore pressures in the soil that 
deteriorate significantly with time. In fact, if ice were present 


at depth 3m. under these thawing conditions, Fig. 5.17 shows that 
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Fig. 5.17 Sample Solution for Pore Pressures at a Soil-Ice interface 
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the pore pressure actually reduces slightly with time as the thaw 
plane progresses through the ice. 

The results from this sample problem suggest that although 
the settlements are obviously considerable, the pore pressure 
conditions in the soil above a thawing ice layer may not be as 


critical as hitherto supposed. 
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CHAPTER VI 
ANALYSIS OF FIELD PROBLEMS 


6.1 Introduction 

In this section attention is directed towards the application 
of some of the methods and the theoretical relationships established 
previously to some practical problems that might be encountered in 
the field. 

One of the most dramatic and indeed topical aspects of 
industrial development in the Arctic is the effect of the operation of 
a hot oil pipeline on the underlying frozen ground. Some important 
effects of a pipeline on its foundation are analysed over the expected 
period of operation for an Arctic location. In the following paragraphs, 
improvements to a thawing foundation using vertical sand drains are 
also considered. 

It would obviously be extremely desirable to obtain many field 
case histories, and anlyse them within the context of the theory of 
thaw-consolidation. It is only possible, however, to present analysis 
of one case study at this time, as there is only one well-documented 
set of data available to the Author at the present. This concerns the 


behaviour of the 48-inch warm pipeline at Inuvik which will be dealt 


with in some detail. 


6.2 Thawing Under a Hot Pipeline 


A practica] study of thawing under a hot oil pipeline is 
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included to illustrate the application of the extension to the theory 
described in section 5.2, where the power law relation between thaw 
depth and time is used. The settlements and excess pore pressures 

under the centerline of a pipeline are analysed, subject to the thaw rate 
given by Lachenbruch (1970). A clay substratum of 65% moisture content 
was considered by Lachenbruch in order to assess some thermal effects 

of a heated pipeline in permafrost. The thaw 'bulb' and the rate of 
thaw under the centre-line are shown in Fig. 6.1. 

Various coefficients of consolidation Cy have been used to 
estimate the performance of the pipe for the first 20 years. The 
material overlying the pipe is assumed to act as a surcharge loading, 
and for simplicity to be comprised of a material of much higher 
permeability than the underlying frozen soil. It therefore does not 
play any active role in the drainage conditions. 

From properties chosen by Lachenbruch, one finds that 2.42 m 
of surcharge material at a dry density of 1.75 t/m°? causes a surcharge 
loading of 4.25 t/m°. The thawing soil at 65% moisture content has a 
submerged unit weight of 0.62 t/m°>. The rate of thaw for the Arctic 


(clay) case is well approximated by 


X mmde5e to: (6.1) 
or B. 2 3256 
and n = 0.3. 


The coefficients of consolidation taken to assess the 


performance of the pipe are 
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Fig. 6.1 Thawing under the Centre-line of a Hot Pipeline 
(After Lachenbruch, 1970) 
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Finally, the closure of the thaw bulb as calculated by Lachenbruch is 
assumed not to affect the one-dimensional drainage under the centre- 
line of the pipeline. 

The results for the normalised excess pore pressures are shown 
plotted with time for different Cy, values in Fig. 6.2, for the centre- 
line of the pipeline. 

If the thaw depth against time curve is approximated by a square 
root of time relationship, and an a value approximately calculated for 
the first two years of thawing, the dotted line in Fig. 6.2 shows the 
excess pore pressures calculated using the linear theory reviewed 
earlier. 

As expected, under these thawing conditions the excess pore 
water pressures are the most severe in the initial few years of thawing. 
For a silty clay with a Cy, of 15° 10°" cm/sec, excess pore pressures 
equal to 50% of the effective overburden loading may be expected even 
after two years of operation. However, it also appears that if no severe 
stability problems occur in the initial years of operation, then 
increased stability may usually be expected as time proceeds, provided 
that soil conditions do not deteriorate significantly with depth. 

The settlement of the pipeline is comprised of the sum of the 
consolidation settlement, and the settlement resulting from the density 
change in the water phase on thawing. The consolidation settlement at 


any time may however be shown to be proportional to the product of 
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Fig. 6.2 Pore Pressures and Settlements under the Centre-Line 
of a Hot Pipeline 


a a 5 te ey 


j 


fom ee cae ee 


. 


X~ 
padi ha wale me) 


a Me) A ae em 


sceViena} “OED es _ 


ls WCRI. MN 
» (POR ete) ara Ovi 


—— 


te 


ray fe ax OES: 


{awoy) t MIT 


acl: ; 2 


Pe a? 


Sah nt! ha ica ‘type ri 


et 


ve eer oT 40: ae 


- OF 


2r3 


S,/S 


consolidation Sf ay is low. At longer times, the velocity of the 


ae arid the depth of thaw X. At early times, the degree of 


thaw front decreases, and S4/ Sia y increases. Hence a more linear 
settlement-time curve is to be expected than that suggested by the 
thaw depth-time relationship. 

Figure 6.2 also shows the transient settlements, Sis under 
the centre-line of the pipe if the coefficient of volume compressibility 
Sc cm*/kg. It is seen that the curvature of these settlement-time 


relationships is considerably less than the curvature of the thaw depth 


vs. time function given in Fig. 6.1. 


6.3 Thaw-Consolidation with Vertical Sand Drains 

In cases where the soil properties in a thawing soil foundation 
are such as to cause the build-up of excess pore pressures and the 
possibility of an unstable foundation, it may be desirable to consider 
certain remedial measures which would improve foundation conditions in 
an economical manner. One of the best-known methods of improvement 
for deep deposits of unfrozen soils is by the use of vertical sand 
drains, (Casagrande and Poulos, 1969; Barron, 1948). In general, the 
presence of a sand drain boundary with zero excess pore pressure serves 
to shorten the drainage path for excess pore fluids, and accelerate 
the rate of consolidation. 

The sand drain problem is usually a two dimensional problem 
in the radial direction r, and the depth x. For unfrozen soils the 
two-dimensional equation may be solved by a method due to Carrillo (1942) 


who showed that the solution to the two-dimensional problem is the 
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product of a pair of one-dimensional solutions. 

Attention is now turned to the problem of accelerating the 
consolidation of a thawing permafrost foundation. Vertical sand 
drains have been installed to improve foundation performance of a 
thawing foundation under dykes constructed on permafrost in Northern 
Manitoba. However, when referring to the design of foundations at the 
Kettle Generating Station in Northern Manitoba, MacPherson et al. (1970) 
Stated that the use of sand drains had limitations because "the extent 
of the drainage required to carry away the liberated water as the 
permafrost thaws is not known". Clearly, it is of considerable interest 
to explore the possibility of solving the required equations, and make 
available a program which might provide a firmer basis for rational 
design. 

It is doubtful whether the solution to this problem is available 
in analytical form, and the method of Carrillo cannot be utilised due 
to the differing nature of the boundary conditions. 

Figure 6.3 shows schematically a typical situation involving 
sand drains, and also shows how the "effective radius" of the drains 
is calculated from the spacing. It shall be assumed that the thaw plane 
remains planar at depth X(t), and that a one-dimensional thawing 
Situation is present, ij.e., 

X(8) =Tay t (6.2) 
As cylindrical symmetry is preserved at all times, the consolidation 
problem is two-dimensional, with vertical and radial flow components, 
but without tangential flow. Accepting that the Terzaghi-Rendulic 


consolidation theory is adequate to describe the flow of water in the 
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thawed soil, the governing consolidation equation is 


atone is Ope ances 4 ida tS aGe 


CLT oodisd 
2 r or (6.3) 
Implicit in this equation are the statements that the 
consolidation coefficient is isotropic, the material properties are 
homogeneous, and that the total stress of any soil element is constant 
with time. 
In designing the numerical procedure for solving this problem, 
a stationary co-ordinate system might be adopted, whereby the finite 
difference equation for the thaw line must be written at a moving 
boundary. In this case, considerable simplifications are introduced 
if a transformation is employed which places the co-ordinate system 
in motion, thus rendering the thaw boundary stationary. 


Setting z = uey the governing equation (6.3) is rewritten as 


or (6.4) 


As the r co-ordinate is orthogonal to the z-direction, the last term 
in the equation is unaffected by this transformation. 


since “4 = 72/7 € eq. (6.4) becomes 
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The derivatives ou and Su bose no special problems in the implicit 
finite difference method adopted, and a set of linear tridiagonal 
equations are still obtained when the finite difference equations are 
written in matrix form. 

The boundary conditions in terms of r and the new depth 
co-ordinate z are summarised as follows. At the interface with the 


sand drain, the soil is free-draining, so 


s 

i} 

m@ 
v 


Ox Zain Us=0; C50 
At the effective radius, from symmetry, no drainage occurs, and 
febs-) 0 < 2 1t0, 2 = 02)  sbese0 
3 3 7 3 
At the surface, a free draining surface is provided, and 


CEN EC gs tt Ve Us=20; en 


At the thaw interface, vertical continuity of the pore water provides 


as before 
ou 
Cy Oz 
At <a Be  Zi eas Al he om + y'X - us : aX? t > 0 (6.5a) 
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Hence the excess pore pressure conditions at the boundaries are 


completely specified. 
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As it is not feasible to obtain an analytical solution to this 
problem, a finite difference numerical method is adopted. 

There are many approaches available for solution to this type 
of two-dimensional problem. Few, however, have been used in solving 
problems with a moving boundary. The Alternating Direction Implicit 
method (ADI) is adopted to solve the differential equation. The method 
is discussed by Peaceman and Rachford (1955) and basically the scheme 
provides a method of numerical integration which does not need excessive 
computation time. For a linear 2-D equation of the heat conduction type, 
it has been shown to be stable and convergent for all values of the 
time interval At. If the subscript i refers to the number of discrete 
intervals in the z direction, and j to the number of intervals in the 
r direction, and k in the time direction, then the method is summarised 


as follows: 


Tog ity ae C sy 
At V aN ed a; SAV, 5) 
+ ey 52( Veer Cy ) (6.6) 
IE a ewe oe CU 7 


where 6° 6 5° and 6, refer to the second and first order 
central difference operators in the r and z directions 


respectively. 


The resulting set of tridiagonal equations is solved for the 
matrix V. ., which is an intermediate set of values at At/e. In this 


first step, the integration is carried out in the r-direction, and 
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using the values Vi the direction of integration is changed, and 
the following set of equations leads to the solution at the (k + 1) 


time step. 
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The equations (6.6) and (6.7) are written in finite aifference 
form according to the Crank-Nicholson scheme. It is still possible to 
obtain a set of linear tridiagonal equations even though the presence 
of the first order derivatives 6, and 6, complicate the resulting 
equations slightly. This is only possible if these derivatives are 
raised to the first power, as is the case here. The details of the 
finite difference equations and the program are given in Appendix A.7. 
The tridiagonal equations are solved by the Gaussian Elimination 
technique used previously in section 5.2. 

The computer program calculates the normalised excess pore 
pressures, and then a vertical integration is carried out on each 
vertical column of nodal points to calculate the vertical settlement at 
each point on the surface. A settlement profile may then be plotted 
as thawing proceeds. 

The program is used to calculate excess pore pressures and 
settlement profiles at the surface for a typical set of conditions. 


The results are shown in Figs. 6.4 and 6.5. 


The excess pore pressure distribution at the effective radius 
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Fig. 6.4 Excess Pore Pressures at 
the Effective Radius 
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Fig. 6.5 Surface Settlement Profiles Around a Sand Drain 
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shows that the excess pore pressure at the thaw line reduces with time, 
and the pore pressure gradient increases. The results suggest that 
the depth of thaw must become significantly larger than the effective 
radius b, before any satisfactory reduction in pore pressure is obtained. 
In the example shown, the drain spacing is 10 feet giving an effective 
radius of 5.25 feet. The excess pore pressure at r = b is reduced from 
66% to 40% after four years, when the depth of thaw X = 20 feet. In 

the same time however, Fig. 6.5 shows that the degree of settlement is 
accelerated from 33% to 77%. So it might be suggested that no 
Significant improvement in foundation conditions need be expected until 
the thaw depth is considerably greater than the sand drain spacing. 

This observation suggests another design alternative. If the 
sand drain spacing is varied with depth in such a way that when the thaw 
depth is small, the drain spacing is reduced, and at greater depths 
the spacing may be greatly increased. This would simply be achieved 
by incorporating two or three different sets of sand drains, each set 
penetrating to a different depth, as shown in Fig. 6.6. The analysis 
of such a foundation might be accomplished in an approximate way by 
changing the effective radius b at certain specified depths in the 
accompanying computer program. 

Certain other practical problems accompany the use of sand 
drains as a remedial measure. If the surface of the sand drain is not 
in contact with continuously thawed soil or water, then freeze-back of 
the surface of the drain might cause drainage impedance and build-up 
of excess pore pressures. However, the sand drain will at least be 


functional for approximately one-half of the year. Also, if large 
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settlements occur in the thawing soil, it is unclear how the drain 
itself would deform, and if it would carry a Significant portion of the 
Structural load on the soil. 

Further investigations are needed, but these preliminary 
considerations suggest that the use of sand drains in a thawing 
foundation might be less effective than hitherto believed, unless the 


design is optimised in a manner such as that shown in Fig. 6.6. 


6.4 Analysis of the Performance of a Hot Oi] Test Pipeline at Inuvik N.W.1 


In order to study the behaviour of permafrost as a foundation 
material for a warm-oil pipeline, Mackenzie Valley Pipe Line Research 
Ltd. installed a 27 m. test section of 61 cm diameter pipe near Inuvik, 
N.W.T. The objectives of the field test were to measure pore water 
pressures, temperatures and settlements in the foundation as thawing 
proceeded. Frozen bulk densities of the foundation materials were 
measured in order to predict the thaw strains underneath the pipe, and 
comparisons could then be made between the predicted and observed values 
for pipe settlement. 

It is proposed to examine the data from the test site in some 
detail, as this constitutes the only completely documented case of 
thawing of a fine grained permafrost soil that is available at the 
present time. 

The data and references used in this study are taken from 
publications by Rowley et al. (1972) and Watson et al. (1973). 

The stratigraphy of the test site is shown in Fig. 6.7. Hot 


oil at 71°C was circulated through the pipe from July 21, 1971, and 
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the results are analysed for the three months following this date. 

About 36 days after thawing had commenced the thaw bulb entered a layer 
of dense gravelly till, and no further significant settlements occurred. 
In the intervening period the thaw bulb moved through a layer of ice 
rich silt, grading from almost pure ice just below the pipe to about 
30% excess ice at the upper boundary of the dense till. 

Settlement gauges and piezometers were installed at the depths 
shown in Fig. 6.7 and the data from this instrumentation is analysed 
here in the context of the thaw-consolidation theory. Figure 6.7 
does not show the horizontal position of piezometers and settlement gauges, 
which were positioned for the most part 0.3 m from the centre-line of the 
pipe. The approximate moisture content of each layer is included. 

The properties of the ice-rich clayey silt layer are given 
by Rowley et al. (1972). The grain sizes are comprised of 53% silt 
and 40% clay, the remainder being sand. Undrained strength values 
obtained from unconfined compression tests and penetrometer readings 
varied from 925 to 1000 p.s.f. at the thaw front, and 1300 to 2000 p.s.f. 


at the top of the silt layer. In situ permeability tests conducted 
~4 4 


on the thawed silt layer gave values from 0.53 x 10 ° to 2.29 x 10° 


cm/s, averaging | x 10°" cm/s. However, data from laboratory tests 
(Watson, 1973) on the same material gave a permeability of 0.17 x 107 
to 0.26 x 107" cm/s. 


Thermal Calculations 


The rate of thaw underneath the centre-line of the pipe was 
observed in the test facility, and this thaw rate is compared here with 


a simple theoretical prediction based on a one-dimensional solution 


of the Stefan type. 
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Figure 6.8 shows the latent heat of the soil plotted with 
depth, and the dotted line indicates the relationship assumed for the 
theoretical prediction, based on moisture content data. 

The expected thaw strains based on frozen bulk density for 
each layer indicate that on thawing, the soil will consolidate in the 
thawed zone to approximately 45% moisture content. From the data of 
Kersten (1949) presented in Fig. 2.5, it is estimated that the thermal 


conductivity of the thawed soil will be given by 


3 


yeas 252) & 108 cal/cm°C sec. (6.8) 


The surface temperature at the pipe base is given as 


From Fig. 6.8 the latent heat may be expressed as 


Ev e/ 0 =O 579" X (6.10) 


where xX is in m. 
From the predicted thaw strains, it is expected that the ice- 


rich clayey silt layer will settle approximately 0.9 m over a depth 


of 2.2lm, therefore 


a hee 
era 

or 
S = 0.403 X (6.11) 


It is not desireable, however, to calculate a depth of thaw 


from the initial position of the pipe, as in this case it would be 
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Fig. 6.8 Estimated Latent Heat with Depth 
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meaningless considering the large movement of the heat source. A depth 


of thaw is defined from the heat source (pipe base) as 


X=X-S (6212) 
and from eq. (6.11), the settlement is re-written as 

5 a ee Kaigll 614 ge (6.13) 
Substituting eq. (6.13) and (6.12) in eq. (6.10) a value for latent 


heat in terms of the new variable X is obtained as 


X= ¥ 46S = ).674-% 
and therefore 


[= RO7% 16:32,.% (6.14) 


Using the simple Stefan formulation, and combining equations (6.14), 
(6.8) and (6.9), a relationship is obtained between the thaw depth 


from the pipe base, X, and tine. 


a Sew 
L(X) 
x= /309t (6.15) 
70 - 16.32 X 


This equation is a cubic in X, but may most easily be solved 


or 


by rearranging in the form 


ihe (22.646 - 5.28 x] x e616) 


where X is in metres, 


and, <. in days. 
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Equation (6.16) is plotted in Fig. 6.9, together with the 
observed values taken from Watson et al. (1973). These values are 


included in Table 6.1. 


TABLE 6.1 PREDICTED AND OBSERVED THAW PENETRATION 


a te SORE NS RL RL oie, 9 oe Ae eee lL GL en SS 


El. pipe El. thaw Date Time time Observed Predicted 

base depth t Wz thaw _ thaw 

(m) (m) (days) days ’“ depth X depth X 
(m) (m) 

1.80 V5 July 22 0 0 0 0 

2.03 2.80 July 28 6 2.45 OL77 Om55 

2.39 Bae) Aug, 5 14 32/4 0.86 0.87 

2.60 B73 Aug 15 24 4.90 Tet3 1.20 

2.0¢ 4.30 Sept 1 4] 6.40 1.48 dis7 

2.98 4.90 Octet 71 8.43 1.92 - 


The comparison between a simple theoretical calculation of 
this form and the observed data in the early stages is excellent, and 
the deviation at later times is possibly due to the two-dimensional thaw- 
ing effects around the pipe. The data could not be continued past 


October Ist as the hot oil supply failed at this time. 


Final Effective Stresses 


The final effective pressures are readily calculated on the 
basis of original moisture content data, assuming 100% saturation. 
The thawed soil above the level of the pipe base is assumed to act as 


a surcharge loading. Although principally comprised of gravel, a 
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Fig. 6.9 Predicted and Observed Thawing at Inuvik 
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small depth of silt and compressed organics which are initially thawed 
also augment the surcharge. The anticipated large strains in the ice- 
rich clayey silt appear to complicate the calculation of overburden 
stresses. However, it may be argued that even though excessive strains 
are taking place above a given point in the thawing foundation, the same 
height of soil solids, and therefore the same final effective stress, 
remains above this point until the thaw plane passes. Thus, the final 
effective stress may be calculated accurately for any point on the 
thaw plane from the initial moisture content data for the frozen 
material. The ultimate effective stresses are calculated in this 
manner from the moisture content data given in Fig. 6.7, and the 
calculations are provided in Table 6.2. The first three layers may 

be assumed to comprise the applied loading Por The remaining layers 
of interest are shown to exhibit an average submerged unit weight of 

y! = 0.434 t/m®. The variation of the final effective stress is 

shown plotted against the original depth from the ground surface in 
Fig; 610. 

TABLE 6.2 FINAL EFFECTIVE STRESSES 


Layer Design- Original Original Unit weight Stress 9 
ation depth (m) height Hm) my Gta inc) YH, t/m 
Gravel ~ O = 1337 1337 122 1.67 tae 
fo) 
Organics ~ 67. = 1452 OL15 0.22 0.03 eh 
Silt clay A Wee roe 0.30 0.22 0.07 
Ice B 1.82 - 1.98 0.16 0 
av. 
I + AS] : - 2, : 
ce + sii Ac 1.98 - 2.74 0.76 0.32 aa Arey 
Sige: fice) 0 2.74 - 3.96 ya 0.51 0.622 
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Fig. 6.10 Final Effective Stresses 
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Fig. 6.11 Settlement Records in the Thawing Foundation 
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Excess Pore Water Pressures 

Attention is now turned to calculating from the recorded data 
the excess pore water pressure at a point when the thaw plane passes 
that point. Watson et al. (1973) present data from eight piezometer 
locations, at depths shown in Fig. 6.7. When the thaw plane passes a 
piezometer tip, the actual head of water is read from the data. The 
position to the water table is monitored continuously from open stand- 
pipes near the surface. The excess head of water may be calculated. 
The excess pore pressure is then normalised by dividing it by the final 
effective stress at that point (See Fig. 6.10). A summary of the data 


1S presented in Table 6.3. 


TABLE 6.3 ANALYSIS OF PIEZOMETER DATA 


2 ; : 
Biezonetan ye bepohel Ue Ctimay | (Ps ey ral ume thy |i el 


(m) ue (t/m2) i 
TP2 2509) 0.3 2.00 Loy Oens 
ei S250 0251 2.42 (aye 0:37 
TP3 2309 0.76 12.99 Be Laas 0.1 
TPS B30 0.42 (We 18.1% 0.31 
ie 4.2/7 0.49 3.00 16.3% Ox? 
TP4 1269 0,25 1.65 Secs 0 
TP4 2.09 0.34 1.98 Vee Ora 
TP4 5250 1.08 2.40 35 .0% Oi35 


The normalised excess pore pressures vary between 15% and 21%, 


excepting two readings which are appreciably higher. 
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Settlement Data 


Two sets of settlement data permit analysis (S2 - S6) and 
(S3 - S10). The gauge $3 was installed at the same level as S2. Other 
sets of gauges were present, but were too far removed from the axis of 
the pipeline to be of concern here. 

The data are best analysed by considering the relative 
movements of a pair of gauges located at different depths in the 
compressible silt layer. Long after the thaw plane has passed the 
lower gauge and entered the dense till, the settlement curve levels out, 
and this value is used for the maximum settlement at that point. This 
value of maximum settlement may then be used to calculate m,, for the 
compressible silt layer. 

In order to obtain values for consolidation settlement, the 
observed readings must be corrected for the following factors. 

(i) "Bowing" or bending of the settlement rods. 

(ii) Thawing of the pure ice layer at the pipe base. 
(iii) The ice/water volume contraction on thaw. 

For (i), the correction is clearly seen from Fig. 6.11. The 
gauge So had moved 0.02 m when the thaw plane passes it, however the 
gauge S6 had moved 0.12 m when the thaw plane passed that gauge. This 
introduces a correction of (0.12 - 0.02) = +0.10 m to all readings 
(S2 - S6) due to bending of the settlement rods caused by horizontal 
movement of thawed material. 

For (ii), the pure ice layer 0.16 m thick and a very small 
depth of soil above the ice layer provide a correction of -0.19 m, 
when making settlement calculations for the silt layer alone. 


For (iii), the ice water contraction results in a volume 
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contraction of 


ay e “frozen ~ “thawed 
1/w 


Ue fe 
et © Frozen 
y 
Now e a ee 
frozen = se thawed 
NG 
2%, G-wQ&G 
(*r Tice : > 
V i. 
1/w y 
NMaris== 4 -G 
Yice : 
ve 4 
(ay) . lice 
et te hs 
ae ace Cn) 


S2 - S6 

The bending correction (i) is 0.1 min this case. The correction 
for the pure ice layer (ii) is -0.19 m, and the correction for the ice- 
water transformation (iii) is -0.0625 Ho: The initial elevations of 
the settlements gauges are 

Si EU .=el.34 

S65 BL =r se00 
therefore the original height of frozen soil is 


SZ = S60" =,2.010.M 


and A) = 0.0625 x 2.21 = 0.138 m 
i/w 


The initial height of compressible thawed soil therefore is 
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Zee) = 05138 =" 0:19 


The coefficient of volume compressibility Mm, is defined as 


5 
ma X 


Mm, = HP. + 0.5y°H) P. a DESy A (6.18) 


Applying the preceding corrections to the maximum settlement 


value in the data, provides 


ma x 
0.632 m. 


Therefore from eq. (6.18) 


0.6319 


My ae NIeSS2F ayy ee Obl Te Tee 


0.86 + 0.1 - 0.19 - 0.138 


= 0.152 .m-/t 


Again consulting the settlement data in Fig. 6.11, the settlement, Sy 


when the thaw plane passed Ses is given by 


S+ =#0,79° 4 09> -".0..19¥=" 02138 


S, = 0.562 m. 
S 

E_ = p06 = 0.889 
max 


Le at the depth S6 


or We = 0.461. 


.. 0.434 x 1.882 
a ean) rf CE 


Using the theoretical curves in Fig. 3.5 for Sf Sha vs. R, the 


required R value is 0.32. 


Combining this value of R with the observed 
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a value, the consolidation coefficient is found to be 


ee 
c= 1.46 x 10° cm/s. [500 ft?/yr.] 


Sa - S10 
Using a similar approach for the relative settlements of 
gauges S3 and S10, the following is obtained 
533 EL = 1.36 
A = 1283 
SiOse EL = 3019 
Sas 0756, SA ag 0.72, and hence the settlement ratio 


Shs 1S" 00/777: 


t’ ~max 


my for the south side of the pipe from this data using the 


corrected height of thawed soil initially as 1.526 m. is 


m, = 0.224 m</ton. 


Calculation of Intermediate Transient Settlement Values 

Using the Mm value obtained for S2 - S6 at maximum settlement, 
it is of interest to predict the transient shape of the settlement 
curve with time. Again, settlement due to the ice/water volume change, 
thawing of the pure ice layer, as well as consolidation settlement 
must be considered. 

Bearing in mind that X is the depth from the heat source 
to the thaw plane, 

X=avt (6.19) 


It is also known that Saty was defined in terms of the initial, 


Levit 008F G2 


0 ag 


i op a ‘ot ar 


otter gnsmelstez shy sored bas (30.0 = 


; a | s st s 
® \ ry * 0 at xem 4° ae i ; " -_ ae 
Le) 


aid pntew arsb atid moet aqiq add to obte dtuoz ot 10% ee 


— ee } 
at’ .m 882.1 a5 gifettint tree oo s fi 
ay 


| var vata 
isi doonstt3o@ snstengaT tei vate 
ednanalt toe ihunitxemr 45 a2 ~ $8 sof Denture Nef 
jnemafdisz sf to -sapite Ses ranayd: ons voRbeig 02 semvaint ¥0 ) 
ie sw fow. 4otaw\Sot orld 03 ok: snaaaidiase ate amit ine v 
suomehaye: noH.abiToenos 28 The 2a ayet spb site: tt 1 


Pi ‘inate i rvs 
a9 102 $5, ‘ond mont yen 


(21.3) 
- yhptsiint silt to emast ab be 


238 


unconsolidated height of soil, thus 


max — m X(P el Cevyisn) (6.20) 
and 
+ 
S = constant (66.21) 
max 
Now S+ Xe Ke oy th (6.22) 


Combining equations (6.20), (6.21) and (6.22) provides the relationship 


>t Dik cay ats 
sua Mm, X(P + 0.5y'X) 
ci + >t 2 
Mm, oe ae Neer: ms Q.5y Xe X= ay bE 
max max 
which gives 
een: >t 
O.5y'm, S Xx” + m,, De § - )X+avt=0 (6.23) 
max max 


For a given time t, this is a quadratic in X, i.e. the depth 


of thaw if no settlement were occurring. 


When the constants for the case in question are introduced 


into eq. (6.23) the following quadratic evolves: 


KCN, (0 ae [162.7 VIAL & (6.24) 


where X is in metres, 
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obtained. 
A table of observed and calculated settlements including the 
necessary corrections for ice layers and the ice/water phase change 


is given in Table 6.4. 


TABLE 6.4 OBSERVED AND PREDICTED SETTLEMENTS 


OBSERVED PREDICTED 
a acne st ee 
days m. m. i/wt+ice 
m. 

eUriyT 0 0 0 0 0 0 0 
23 ae Als: 1.42 0.427 Onli2z 0.109 0.330 
Zo a ed 2.00 0.605 ORSAY 0.157 0.392 
27 6 0.23 2.45 0.744 Oe) 0.197 0.442 
29 8 0132 2203 0.863 0.260 0.231 0.485 
3] 10 0.39 3.16 0.966 0.295 0.262 02523 
2/8/71 12 0.44 3.46 1.060 0.326 0.290 0.563 
a 14 0.53) 3.74 1.149 OF357, 0.318 0.593 
6 16 0.56 4.00 lee ol 0.336 0.343 0.624 
8 18 0.60 4.24 1.308 0.414 0.368 0.654 
10 20 0.65 4.47 1.381 0.440 0.391 0.682 
i 22 0.70 4.69 1.452 0.466 0.414 Onn 
14 24 0.74 4.90 lo 0.491 0.437 0739 
16 26 0 ,2.5 sam OLG10 1.584 0.516 0.458 0.764 
18 28 0.80 5129 1.645 0.539 0.479 0.790 
20 30 0.83 5.48 UB ZACH 0.563 0.500 0.816 
a $Y4 028551" 5:66 PhO 0.585 0.520 0.840 
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Table 6.4 Continued/... 


ee eee sey) Wee a LES oR ees) 2 vias ale ee 


OBSERVED PREDICTED 
Date c S vt rn = 
days m. Wer ink max t : ze ; 
days m. m. i/w + ice 
m. 
24 34 OF37 5.03 1.821 0.607 0.540 0.864 
26 36 0.89 6.00 1.876 0.629 0.359 0.887 


24/9/71 64 0.96 8.00 ~ - é e 


Fig. 6.12 presents the predicted and observed data for (S2 - S6). 
The curves must agree at the "end thaw" time, as m, was calculated from 
the observed data in any case. The real value of the plot is the 
comparison between transient settlement behaviour prior to complete 
thawing. The predicted values agree well with the observed settlement 
curve at later times. In the earlier stages, however, the theory 
considerably overestimates the settlement rate. This may be because in 
the theoretical treatment, the 0.16 m ice layer was presumed to thaw 
almost instantaneously, as indeed probably happened. However, the 
field settlement plot does not reflect this sudden settlement to the 
same extent, and consequently lags behind the predicted values. 

The reasons for this are not clear, but it may possibly be due 
to non-uniform thickness of the ice layer, giving an arching effect 
in the overlying soii. Horizontal movements of material slipping in 


towards the pipeline may also have delayed the vertical settlements. 


Summary of Results 


The Bip) > ase values obtained from the two sets of settlement 
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data enabled two calculations of the required R value in the field to 
be made. The required coefficient on consolidation was then calculated. 
As the Mm, value for each location was known exactly from the total 
settlement data, the permeability k required to give the degree of 
consolidation could be calculated. These two values of the required k 
are included in Table 6.5, together with the data and similar calculations 
for the eight piezometer locations considered. 

In general, the theory of thaw-consolidation indicates that 
the excess pore pressure data and the settlement data are entirely 
consistent with each other. This is so because the required 
permeability values to obtain the observed settlement ratios and the 
observed excess pore pressures lie within the same range of values. 
The theory would indicate that the following narrow range of soil 


properties were exhibited in the clayey silt layer on thawing: 


Pa) y-ee nee. 10°* om*/s (400 - 600 ft¢/yr) 


m 1.7 - 2.2 cm°/Kg 


i oe OL Ao se ar 


cm/s 

It now remains to compare the required permeability in the field 
with those that were actually measured for this material in field 
and laboratory tests. 

Rowley et al. (1972) have given values for permeability obtained 


4 4 cm/s. These values seem 


aed. field testeof C.53°x 10, to 2679 % 10" 

high for a soil of this type, and it was suggested that the macro 

structure of the thawed permafrost was responsible for the large values. 
An open stand pipe with a 9 inch by 1 1/2 inch cylindrical 


tip was used in the test. The head of water used in the stand pipe 
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was approximately 3.0 m (10 feet) above the ground water table. The 
tip was brought close to the thaw plane at different times, and the 
test carried out under a constant head. This head of water represents 
a water pressure of 0.3 Kg/cm* in excess of the hydrostatic ground 
water condition. Referring again to the diagram of final effective 
stresses with depth in Fig. 6.10, it is seen that at no point in the 
clayey-silt substratum does the effective overburden pressure approach 
a pressure of 0.3 Kg/em?. This means in effect that the water pressure 
exceeded the maximum effective overburden stress, causing hydraulic 
‘jacking' of the substratum, and the permeability results so obtained 
must be considered suspect. Field permeability testing in such 
Situations must be carried out with great care, giving due consideration 
to the range of effective stresses involved. 

Laboratory data on the permeability of the clayey-silt layer 
have been provided by Watson (1973). The soil samples were thawed 
under a small stress, and the permeability calculated from constant 
head tests. The void ratio and permeability were determined for 
Successive increments of loading. The data from two such tests have 
been provided, and the results are summarised in Table 6.6 


TABLE 6.6 LABORATORY DATA FOR INUVIK CLAYEY SILT 


Load o' (Kg/cm*) Void Ratio e Permeability k (cm/s) 
ROR TTR rene Ee ie aes oe ee 
0 Sed - 
-4 
0.186 We Oe24 2x) 1.0 
0.293 1.16 ipatelce 10 
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0.880 1.02 0.043 x 10° 
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TABLE 6.6 continued... 
een eee re A eee en ET CRO RY pole Wil a 


Load o' (Kg/cm*) Void Ratio 3 Permeability k (cm/s) 
Ren AT an te ae a foe ae es 
0 3.01 : 
0.293 1.22 9.108 x 10°" 
0.585 1.09 0.052 x 1077 
0.880 1.02 0.027 x 1074 


Figure 6.13 shows the void ratio plotted against effective 
Stress, and the permeability, k. Both relationships are linear on a 
semi-logarithmic plot. Using these relationships, the permeability 
may be plotted directly with effective stress, to obtain the graph 
shown in Fig. 6.14. 

The graph also shows the narrow range of effective stress 
experienced by the clayey silt layer shown in Fig. 6.10. Taking the 
widest range of permeabilities possible, the permeability of the 
laboratory samples in the stress range of interest lies between 


" : cm/s. Although this range of values is 


OFl7 2410) -« and 0.260% 105 
determined from specimens of the foundation soil tested in the laboratory, 
and it is true that the macro-structure in the field may not be 
adequately represented, the laboratory values are considered a great 
deal more representative of the soil behaviour than the results of 
the field tests, which must be regarded as too high, for the reasons 
cited above. 

Table 6.5 indicates the use of the observed pore pressures 


and settlements to obtain the R value in the field that is required. 
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for consistency with the theory. As the a and Mm) values may be 
accurately calculated from the observed performance of the foundation, 
the remaining unknown quantity in the R value is the permeability, k. 
Table 6.5 also shows the calculation of the permeability value required 
in the field to verify the predictions of the thaw consolidation 
theory. There are ten such "required" permeability values at different 
depths, and these are plotted against depth in Fig. 6.14. Table 6.5 
also shows the average required permeability is calculated to be 


0.25 x-fo"* 


cm/s. 
The narrow range of permeability values obtained from 
laboratory testing is also indicated in Fig. 6.14. These values have 


4 cm/s. 


an average of 0.22 x 10° 
The range of permeability values required to validate the one- 
dimensional theory of thaw-consolidation presented in this thesis 


4 4 cm/s. 


is 0.1 to 0.4 x 10 ° cm/s, with an average of 0.25 x 10° 
Laboratory testing on thawed samples at the effective stress range 
experienced in the field gave a range of permeability values of 0.17 


4 cise tee 1 Ss 


to 0.26 x 107" cm/s, with an average of 0.22 x 10° 
therefore concluded that the agreement between the theoretical 
predictions presented here and the observed behaviour of the 
foundation is excellent, particularly when the extreme variability of 
a parameter such as the permeability is recognised. 

It is suggested that if the quantities involved in the 
calculation of the R value had been estimated in advance of the 


operation of the test pipe line, and the one-dimensional thaw- 


consolidation theory used to predict the pore pressures and settlements 
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under the structure, remarkable agreement between theory and 
observation would have been obtained. Both theory and observation 
demonstrate that although the total amount of settlement might be 
considered excessive in some cases, the soil properties were such as to 
ensure a reasonably stable foundation while thawing was taking place. 
Approximately eighty percent of the maximum available effective stress was 
attained by the soil immediately on thawing, and between eighty and 
ninety percent of the settlement took place at the same time. 

The analysis of the field data presented here assumes the existence 
of a homogeneous soil profile, and the maintenance of one-dimensional 
thawing conditions. Clearly, neither of these assumptions are completely 
realised within the thaw bulb. The soil profile exhibits a decreasing 
ice content with depth, and therefore is non-homogeneous in this respect. 
However, calculations of My and therefore of Cc, are based on an average 
total settlement throughout the layers of interest. In this way the 
differing soil properties are averaged, so that the theory for homogeneous 
soils may be used for analysis. The assumption of one-dimensional 
compression and drainage conditions is assured under the centre-line 
of the pipe. Examination of the shape of the thaw line and the surface 
settlement profiles given by Watson et al (1973) assures that at 

Although this test pipe line loop constitutes the only 
complete case history available at the present, it is extremely well 
documented and reported. The comparisons between theoretical and 
0.3 m from the pipe axis (where the instrumentation was installed) 
the assumption of one-dimensionality is still reasonable. 

Although this test pipe line loop constitutes the only 
complete case history available at the present, it is extremely well 
documented and reported. The comparisons between theoretical and 


observed behaviour considerably increases the confidence in the 
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CHAPTER VII 
CONCLUDING REMARKS 


In this thesis, the equations governing the consolidation of 
thawing soils have been presented in a form sufficiently general for many 
applications. In situations where a uniform frozen soil is thawing 
due to the application of a constant surface temperature, a useful 
solution in closed form becomes available. In such instances, the 
thaw-consolidation ratio, R, enjoys a privileged position, as it 
controls to a large extent the maintenance of pore water pressures 
and the degree of settlement. With careful determination of tne 
coefficient of consolidation, it should be possible to provide a 
fairly reliable forecast of the foundation conditions. 

It has been demonstrated how numerical solutions may be used 
to account for features which would have to be excluded or over- 
simplified using an analytical solution. 

The concept of the residual stress was introduced to describe 
the conditions in the thawed soil if no drainage were permitted. 

This parameter has been successfully measured for reconstituted and 
natural samples of thawed soil. Its influence on the build-up of 
pore pressures and the subsequent deformations in the thawing soil 

is clearly important, and some implications are cited for cases where 
the residual stress is large compared with the overburden loading. 


The predictive power of the theory has been shown to be 
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Substantial by carrying out carefully controlled laboratory tests, 
and by analysing in detail the performance of a thawing foundation 
under a hot oil pipeline. 

While the theory in its present form lends itself admirably 
to many applications, the present work indicates that further research 
in certain areas is of high priority. It has been demonstrated quite 
conclusively that at the present state of understanding, further 
research into the thermal properties of soils must be considered 
subordinate to the requirement for a more exact knowledge of 
geotechnical properties. The coefficient of consolidation is the 
most important single parameter in this respect, and it in turn 
contains the permeability of thawed soil. When applying the theory 
to field problems, the correct in-situ permeability becomes a most 
important auaneaies indeed. In large scale projects, the in-situ 
measurement of permeability may be warranted, and it is of paramount 
importance that the permeability be determined at the correct 
effective stress level, as the strong dependence of permeability on 
effective stress or void ratio is well known. 

The preliminary measurements of residual stress described in 
this work have suggested that a significant increase in this quantity 
with depth may be expected, where the void ratio decreases with 
depth. Further testing of samples from natural permafrost profiles 
appears to be of the highest priority, as it may provide great insight 
into the natural structure of permafrost, and assist in the more 
economic design of warm structures in contact with frozen ground. 


The analyses performed here have been limited to establishing 
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the physics of thawing soils in a one-dimensional geometrical 
configuration. However, many field problems are obviously multi- 
dimensional in nature. Although no serious conceptual difficulty 
Should arise in their treatment, multi-dimensional problems might 
very usefully be formulated and solved to augment the analytical 
power of the design engineer. 

Finally, the extreme shortage of well-documented field case 
records is obvious, and any research carried out in this area will 
certainly be well-received, and will greatly enhance the present 
state of knowledge concerning the realistic behaviour of thawing 
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APPENDIX A 


DETAILS OF FINITE DIFFERENCES PROCEDURES AND PROGRAM LISTINGS 


A.1 Finite Difference Procedure for the One Dimensional Thawing of Soils 


This procedure is very similar to that described by Ho, Harr 
and Leonards (1970). The finite difference scheme described here is 
designed to investigate some of the effects of temperature and phase 
properties on the rate of thaw and the temperature distribution in a 
melting soil, and is not intended as a user-oriented computer program. 

The governing differential equation of one-dimensional 


conductive heat transfer is from Section 2.3 


Core (Ry fen 30 | _ L w Yq d Wh 
A) c, (8) ox ox c,(6) dt Caray) 


Assuming a constant discrete interval in the x-direction, and writing 
the governing equation in finite difference form according to the 


Simple explicit method, we obtain 


=~ (Aj+iB snd)eGe att tBedtis 


nae Tiss lel be 


PNM alee} ania 


-Lw 1G AW, /C,(8; 3) (A.2) 


where i andj are the finite difference subscripts in 


the x and t directions respectively, 
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(A.3) 
(A.4) 
and where c,(8, ;) and k (0, 3) represent the values of those 
functions at the temperature 6 = 0. j 
For a stable solution the inequality 
Alas Baxi] (A.5) 


must be satisfied, and this limitation on the time step At is the 
chief drawback to this numerical method. However in this instance, 
this restriction is easily counterbalanced by the ease of 
programming this simple method. 

The quantity AW is the fraction of the total moisture 
content that changes phase in the time interval At. 

Thus 


AW = W 08; 


ina f 1) 5. W C8; 5) (A.6) 


J 


The quantity Wy (85 aan 1) is unknown at the time level j, 
hence the finite difference equation is essentially implicit in form. 
For arbitrary functions of W (8) we must solve each finite difference 
equation by the Newton-Raphson Iteration method. This method is 
extremely efficient, and converges rapidly to give the solution 85 Ao ae 


The arbitrary functions Co(8). k(@) and W (8) are included 
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in the program as function subprograms, and they may be varied by 


changing the required cards in each function subprogram. 


A.2 Program Listing for One Dimensional Heat Conduction 
Purpose:- To calculate the temperature distributions and depth 


of thaw in a saturated soil. 


Description:- The program uses a simple explicit finite difference 
procedure to calculate the conductive heat transfer in a melting soil. 
The latent heat of the soil is accounted for in the manner described 


by Ho, Harr and Leonards (1970). 


Usage:- The data is input on a single data card FORTRAN format (9F8.0) 
in the order 
k 


Say Pa AAS HS th rat 


TR ae g 
where ky and ke are the conductivities of the fully thawed and 
fully frozen soil respectively, 

P, Q and R define the unfrozen moisture content/temperature 
relationship by the equation W, = (P + exp (Qo+R) }/100 
and 6 is the temperature, 

H is the height of soil considered, 


ie and Tg are the surface and ground temperatures, 


and w is the water content expressed as a fraction. 


Units:- The units used are centimeters, grams, seconds and 


degrees Centigrade. 


Qutput:- The program will print out the time in seconds, hours and 
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1/2 


hours and the temperature profile at every time step. 


Example of Data Input:- 
k= 0.00245 cal/em.s.°C 
ke = 0.005 cal/cm.s.°C 


P atade 0 

9 Miog 

Ri Sewer, 6 

H = 500 cm. 

Tee werdc 

Tj tae -3°C 

atonal OMG; (See overleaf for program listing. ) 


A.3 Finite Difference Procedure for Thaw Consolidation with Arbitrary 
Movement of the Thaw Plane 

The finite difference grid is fixed in space in the x-direction, 
and the moving boundary is denoted by node 'n'. The previous (n-1) nodes 
are equally spaced at Ax. For any arbitrary specification of the thaw 
plane 

i Paenex sc) (A.7) 
the integer number 'n' is calculated by 


7 Renee i (A.8) 


and this indicates the position of the thaw plane in the finite difference 
grid. This value is stored in the computer program in an integer 
location, and so any fraction of a grid spacing is truncated. We 


denote the fraction of the grid spacing between nodes (n-1) and n by V, 
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73 


42 


62 
63 


45 


47 


i 


15 
16 


REAL U(200) eUNEXT(200) 

COMMON PeQeRe TCTe TCE oWe GAMO 
READ(Se2) TCT TCE sPeQeReH, TS, TI) W 
WRITE (6083) TCT»TCF oH 


FORMAT(® *.//7% THERMAL CONDUCTIVITY IN THAWED SOIL IS*.,F10.5/ 
* THERMAL CONDUCTIVITY IN THE FROZEN SOIL I1S*%.F1005/ 
THE HEIGHT OF SOIL CONSIDERED [S® oF 10.3) 


FURMAT(3£ 8.0) 
WRITE(6073)PeQsR 


FORMAT (*0°% *P=* eFS8e2% ° Q=*sF Be2, ° R=* ,F8e2) 


N=40 
N=30 

OX=H/N 

GAMD= 2e7/(1e4W2e7) 

CCF=TCF 

IFC TCTeGT.TCF) CCF=TCT 
OT=GAMD#OX#DX *( 0024005 #W)/(2.2eCCF) 
OT0=062*0T 

OTF =O0T 

WRITE (6042) DT .DX 

FORMAT(*°O%.* THE TIME STEP IS*.F10030° 


OO 3 I=1.N 
UCTI)=TI 

OUT=0e4 

TOUT=0e4% 

TIME=06 
TOL=0.00005 
UC(N#1 )=U(N-1) 
TF(UC1))61261.62 
OT=O0TO 

GO TO 63 

DT=O1F 
BET=OT/(OX*DX) 
00 S I=1.N 
T=U(1) 

TP=U(I+1) 
IF(T-EQe1) GO TO 6 
TM=U(I-1) 

GO TO 7 

TM=TS 


A=2e*BETS(CO(TI¥CLeSTC(TP) + Le STCCT))) 


IF(1eGTel ) GO TO 44 
B=26¥* BETH#IC(TM)/CO(T) 
GO TO 45 


B=2e*BET/S(CO(TI&CLe4TC(T) + LeZTC(TM))) 


AB=A+B6 
IF CABGeLT. 10) GO TO 46 
WRITE (6047)1 


FORMAT (%0% » & THE® »13e*TH NODE IS UNSTABLE » 


BE HAL VEDeece® ) 
OT=O0T/1-59 
BET=BET/1-S 
GO TO 4 
E=79e6%W*GAMD/SCO(T ) 
D=A&TPt (1 e—-A-B) HT +B*¥TMFEER WUT) 


DX=* eF1063) 


CCMMENCE THE ITERATICN TO FIND U(ToJdtldeoe 
LET THE INITIAL APPROXIMATION GE UCT oJ) cece 


KOUNT=0 

SOL=T 

STORE=SOL 

F=SOL-O+tE* WUC SOL) 

IF (SOL e*LT 200) fe} Ae) US 

IF (TLT O02) GO TO 50 
FDASH=le 

GO TO 16 

FOASH=1~. + E* (WUC SOL) -WUCT))/ZCSOL-T ) 
GO TO 16 
FOASH=1Le tE# O8EXP(QFSOLERIS1IO0 © 
SOL=SOL-F/FDASH 
DIFF=ABS(SOL—STORE) 
KOUNT=KOUNT #1 
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IF (KOUNTeLYT e20) GO T0 9 

IFCT.GT.-0062) GO TO 105 

SOL=0-5*(SOL*#STORE) 

WRITE(66106)1 

FORMAT(® ®©,® THE SOLUTION FAILED TO CONVERGE AT NODE*,15) 
GO TO 12 

SOL=0.0001 

WRITE(6.10)1 

FORMAT(® 5% THE *%ef4se*TH NOOE HAS COMPLETELY CHANGED PHASE®) 
GO TO 11 

IF(OLIFFeGTeTOL) GO TO 8 

UNEXT(1)=SOL 

CONTINUE 

TIME=TIME+OT 

HOURS=TIME/ 3600. 

SQHR=SOQRT (HOURS) 

OO 12 1=1.N 

UCT) =SUNEXT(T) 

WRITE (6Gc13) TIME SHCURS eSQHRe(UCI) eo 1T=1.2N) 

FORMAT( ® eo *TIME=*® oF l2ele® HOURS=* oF 100e20% ROOT TIME=* eF10.3/ 
(10F10.4)) 

IFCUCN) -LT.-0.1) GO TO 4 

GO TO } 

END 


REAL FUNCTION TC(T) 
CCMMON PeQeRe TCT. TCE &e GAMD 


THIS FUNCTION EVALUATES TKE CONDUCTIVITY AS A FUNCTION OF 
TEMPERATURE eoewe 

SEKKEKKE KEKE HK EK SEKER PEEKS KSSEE KKEKKERERE 

TFC ToL TeO-) GO TO 1 

wU=1e 

GO TO 2 

WU=(PHEXP(Q*T4#tR) )/1006 

TC=TCFewU*x( TCT-TCF) 

RETURN 

ENO 


REAL FUNCTION CO(T) 
CCMMON PeQeReo TCTs TCF e We GAMD 


C—-----THIS FUNCTION EVALUATES THE APPARENT VOLUMETRIC SPECIFIC HEAT 


2 


OF THE SOILeece 

KEKE SEKRKEKKSSEKKEEK KK EERKESESEE EE KKKKKESEKERE 
TE(T LT e002) Go TU 1 

wU=1. 

GO TO 2 

WU=(PtEXP(Q*T#R)IIS100~ 
CO=GAMD¥( 062400 S#we( 1 ot WU) ) 

RETURN 

END 


REAL FUNCTION WU(T) 


ESI I THIS FUNCTIGN EVALUATES THE UNFROZEN MOISTURE CONTENT OF THE SOTL 


(CSS AT THE TEMPERATURE *%4T%*% ececcce 
C----- KRETRKEEKKEKEREK SR HREKRREKE KE SRESKEK SEE KEREKKEKKKEERKESE E 


1 
2 


CCMMON PeQeRe TCT eo TCh oe We GAMD 
IFC TeLTeOe) GOTO Ry 

WU=1 6 

GO ¥O 2 

WU= (P#HEXP(Q#TER)I 71006 
RETURN 

END 
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and it is calculated from the relation 


v= 4 - (n-1) (A.9) 


Thus there are (n-1) nodes spaced equally at Ax and one node with a 
spacing of V(Ax). 

The finite difference equations for the first (n-2) nodes were 
written as usual according to the Crank-Nicholson (1947) scheme for 


the unknown values of the pore pressure at the (j + 1) time interval: 


a CR eect LEAR ae ot bday 
= 8 Uio1, 3 + 2(1-8) Us + 8 User Ge i= Vf, = n=2) (Ant0) 
Cy At 
where 8g = 7 
(Ax) 


At is the time increment, 
i is the nodal subscript in the x-direction, 


and j is the subscript in the time direction. 


The pore pressure at the surface is always taken to be zero, 


and so 


Yo,5 7 Yositl * oe 


For the (n-1) node, the unequal spacing of the nodes requires 


a different equation. 


Let 
Vie] = fractional node spacing at the jt+l time level, 
and V = node spacing at the i level. 


These values of V are calculated from eq. (A.9) using 
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X(t 41) and X(t.) respectively. 
Using the Crank Nicholson scheme 
ye 4 


ou _ u 0 U 
7 1/2 Cy, Ae + aw 
j Soa 


Q 


-) 


we obtain after application of the usual central difference operators 


-B V. 
jt] Pi: shee 
BE Vis 
ou J B 
ot race rice eer omelet Omen Unt 
ie VeD Tidy nay ee nd (A.12) 


For the derivation of the equation at the n th node we must 
assume for consistency with the previous equations the existence of a 
fictional node (n+1) at a distance V(Ax) below the boundary node n. 

At the moving boundary on node n the following equations must 


be satisfied at the two time levels j and (jt+1). 


AL tet, tay, Ste cy a es: 
and Po + y'X - u = Sv Be an 
dt 
Writing eq. (A.13) in finite difference form 
“Bigg Uner ger + 248544) Yn jet 7 Bje1 Yn-1, 541 
Or ie ome Lead Ga Sele (A.15) 
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CL 
where 8... = Sa ea 
Vay (Ax) 
: Cy At 
an B. = Marie 
J ve (Ax) 


Writing (A.14) in finite difference form 


Vas) 
: : j dx 
Unt sj 2(P. Une) pect j Bor en bee (A.16) 
and Une], gt] 
3PM Sot) Pa delay | On-thae (A.17) 


where P and Pang are the values of iP + y'X) at the j 


and (j+1l) time levels, 


and dX are the known velocities of the thaw 


dX} and dX 
dt 


t jt 


J 
fronu ac t and teat 
These values may be obtained by differentiation of the function X(t). 
Equations (A.16) and (A.17) may now be substituted in eq. 
(A.15) to eliminate the unknown values at the fictional node (n+1), 


and finally the equation for the node n is written as 


(1+ Biyy + Bia yCjag) Yn ged 7 8547 Unt, 341 
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V (Ax) 
where C — * 
V J+] (A.19) 
V. (Ax) 
and C. = Ji # 
J « j 


The values for P, C,; 8 and dX/dt are calculated for each new time step, 
and used in eq. (A.19). 

The equations may now be written in tridiagonal matrix form, 
and solved by a very efficient algarithm derived from the Gaussian 


Elimination technique for solving simultaneous linear equations. 


gee She dP ed has th 
Co ay bo Xo do 
Ci a b, x d. 

ay Dy Xn qn (A.20) 


Briefly, if the coefficients of the diagonal, super-diagonal, the sub- 
diagonal, the right hand sides and the solution vector are denoted by 
ass b.5 C.5 d. and Xs respectively, then the algorithm for solution 
proceeds as follows. 


1. Change into upper triangular form using 
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po3) 
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2. Back Substitute using 
X = 0 


x, = (di - be xX5,))/ass i = nm, (n-1), ----,1 


The solution vector is now available in Xs 5 and this will 
generally form the solution at the (j+1) time step, ee | 
A.4 Program Listing for Thaw-Consolidation With Arbitrary Movement 


of the Thaw Plane 


Purpose:- To determine the excess pore pressures and degree of 
consolidation in a thawing soil, when the movement of the thaw plane 


is specified by an arbitrary function. 


Description:- The program uses a Crank-Nicholson finite difference 
procedure to solve the equation of consolidation subject to the 
boundary condition at the thaw line. It uses a fixed grid, and the 
boundary moves across the grid with time. The listing shows a thaw 
boundary movement described by 

Xo bee 

where 8B and n are entered as data. If one wishes to specify an 
alternative function, then the two cards specifying X(t) and oe 
may be changed in subroutine THAW. The card specifying TO on line 24 


of the main program must also be changed to be consistent with the 


initial thaw depth XO. 
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Usage:- The program requires data on two cards as follows. The first 
cards in FORTRAN format (7F10.0) accepts the following data in order, 

Bs. Ms tes Cys AX, Pos yi 
where B andn are parameters describing the movement of the 

thaw plane, 

t. is the termination time for the program, 

Ax is the grid spacing in the x direction, 

Res is the surface applied loading, 

and y' is the submerged unit weight of the soil. 

The second card accepts one parameter in FORTRAN format (12) denoted by 
NG. This parameter changes the size of the time increment used, and 
usually lies between 20 and 200. It may be varied to determine the 
effect of the size of the time step. If it is specified too small, the 
time increments become very larae and instability of the numerical 


procedure may result. 


Units:- The units used are lbs, feet and years;or tons, metres and years. 


Output:- The program prints out the thaw depth, time, time!/, time 


factor, degree of settlement and the normalised excess pore pressure 
distribution, and this information is printed out two hundred times 


between t = O and t = te. 


Example of Data Input: - 


Beas 7 feet/year?** 
ie 4 a0 3 

te = 20 years 

c= 100 ft@/year 


$21tt aT .ewollo® 26 2b1s9 ots aia eovtoper T =39 
etebvo ni ated dntwot oF ant etasoos (0. oF at), 3 80 oan : , 

NY ag 9 ath es | mie a 

eit to Semeyon srt entdtraeeb mae : 

consiq want ee au 

die tno git tat omit noi tan tes ait at 4 

efoh32e1 Fb. %ad3 nt entoege bine ‘od at 

ep rpsof botlags agin ats: fs 

.f toe odd to. grip tew J tnw bagtendue st - a — 

vd batensab (ST). Tpnrv0? PAE TRS tek 18d mR ano sesh agead gett aa 

bas <beaw tnameiont ants oft to este art 2epesi> swig th 

odt snimvetsb ot betyey o¢ yeu $7 00S bre 0s nsondod 291 | 

.Hfsme cod Bertiosge af FP 41 gate smtt ent Yo oste ae 10; 

featysmun srt to. sitll bas aieed rey and 2 

dain, 


Bs 
4 iy rae 


X6 


e 


1H9y bits Baytom ,2nod 1OL2189Y bas toot sia i bind 


sweestg 9169 2a99K9 eis, out bos tants to estpsb rota? 
eomty beri Gws INO botntig en noiteana'tni’ aiis bs l oe ae; a sa 
| ee 


er leant blag Hees 
dat at i a 


| en ” 


SS - J. —_ 


269 


Ax = 0.5 feet 


P= 870 Ibs/ft® 
Sieyimarmnncr fecty it 
NG = 50 (See overleaf for program listing.) 


A.5 Numerical Solution for Thaw-Consolidation in a Two-Layer Profile 
We assume the solution for thawing and consolidation in the 
upper layer to be defined accurately by the analytical solution given 
in Chapter 3. We require a solution for the excess pore pressures and 
settlements when the thaw plane enters the underlying layer at depth H. 
Letting the finite difference interval in the x-direction be AXy> then 
there will be 'm' finite difference nodes in the upper layer defined by 
m = H/Ax, (A.21) 
The movement of the thaw plane in the lower layer is given by 
the relationship defined in Section 2.9. In the same manner as in 
Appendix A.3, the position of the thaw plane in the finite difference 
grid is given by 
X=H+ (n -m- 1) AX, + V(Ax,) (A.22) 


Writing the equations for the top layer we have 
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C----—PROGRAM TO CALCULATE THE PCRE PRESSURES IN A THAWING SOIL. 
c----- THE PROGRAM USES THE CRANK-=NICHGLSON FINITE DIFFERENCE METHODes 
C-----1T 1S A FIXED GRIDs MCVING BOUNDARY SITUATION, AND THE SOIL IS 
C-----ASSUMED FOMOGENEOUSese. ONE DIMENSIONAL FLOW [S PERMITTEDe AND 
C-----THE PROGRAM CAN INCORPERATE ANY SPECIFIED RATE OF THAWeeee 
C-----— EKER SERHKEKRHEKEKEKEKKKCRE KEKE KEKE RES EKER EKRHSEKEKRSEEKEKKEK HERRERO EKE KEKE 
REAL U(S500) 6A(500).8(500) 2€(500) .D( 500) 
CCMMON DXsPO 
C------ ENTER THE PARAMETERS NECESSARY FOR SOLUTION 
99 READ (56101) CBeCNo TF oC Ve DX 0 PO. GAMD 
101 FORMAT(7F1000) 
WRITE (66102) Cte CNo TF eCVeDX sP0 eGAMD 
102 FORMAT(*1%,°® DATA INPUT: -"78 THE RATE OF THAW IS GIVEN BY 
L2°/20Xe"X(T) =" sF6e20% eTIMEt* sF563/°THE SOLUTION WILL BE TERMINAT 
2EC AT TIME="*,Fi0e3/*THE CUEFFICIENT OF CONSOLIDATICN I1S*eF19.5/ 
3*THE GRID SPACING IN THE X-DIRECTION 1S*of7e3/*THE APPLIED LOAD AT 
& THE SURFACE IS*eF10e4/*THE SUBMERGED UNIT WEIGHT OF THE SOIL IS*s 
5F 8.3) 
REAC(5e31)NG 
31 FORMAT(I2) 
TOUT=TF/2006 
TOUTP=TOUT 
XO=2.1#DX 
TO=(X0/CB) *€#(16/CN) 
TIME=TO 
UC1)=062*P( PO »GAMD, x0) 
u(2)=U01) 
DT=TO/20. 
X=X0 
DERIV=X/TIME 
CALL POSITCXsDXeNeVeNP) 
N1=NP 
Vi=v 
N=N141 
UCN)=U(1) 
CALL GUTPUT(UsX eN1l oe V1 sGAMD oT IME) 
25 TIME=TIMEFOT 
IFC TIME sGT.TF) Go TO 99 
BT=CV*DT/(DX*DX) 


c----- STORE xX AND DX/DT AT THE JTH TIME STEP eccce 
X J=X 
DER J=DERIV 
c----- CALCULATE X AND DX/DT AT THE (J#1)TH TIME STEP ecece 
CALL THAW(X eDERIVeTIMEsCBeCN) 
KIT=* 


OERJI=OERIV 
CALL POSIT(XeCXsN2eV2eNP) 
N2=NP 
IREM=0 
IFCNZcEGeN1 ) GO TO 1 
TREM=Ne2-Nl 
V2=V2+tIREM 
MAR K=1 
1 O19 )=26*(1.2-BT) *€U(1 FET F#UC2)) 
N=Ni 41 
Ni 2=N-2 
IFC AL2-LEo1) Go TO 4 
00 2 f=2sNL2 
2 OCT )=BT*¥(U(CT4+1)+4+UCI-1)) #2eo%(1 o-BTI*UCI) 
4 D(N—-1L)J=HBTFEVIFUCN-2)7C1LetV1) #(10e-B8T)*UC(N-1) #¢€ BYXUIN)/( 1 -#V1) 
OO 3 [=1-eNtL2 
ACT )=2¢*C1 -+8T) 


ECI)=-8T 
3 c(1)=-BT 
C(1)=06 


ACN-1)=164UT 
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CCN—1)=-BT¥V2/(1.¢4¢V2) 
B(N-1)=-B8T/(12+V2) 
C~-—--—-CALCULATE THE CUEFFICIENTS FOR THE FINITE OIFFERENCE EQUATION 
Cm-——---AT THE BOUNCARY( NTH) NODE ecoccse 

BI=CVFEDT/(V1I¥V1L*OX*DX ) 
BJII=CV#*DOT/(V2*V2RDX*¥DX) 
PJ=P(PO0eGAMDo XU) 
PJII=P(PO0eGAMDeXJ1) 
CJ=V1I*OX*FDERI/SCV 
CIULI=AV2COX*DERII/SCV 
A(N)=10¢4B851 #(1Le#CJl) 
C(N)I=—-By51 
BC(NI=O06 
D(N)J=(1.-BJ-BI*€CI ) RUC N) CBJ B8U(N-1) #¢B4U1 #CUL*¥PIL +8659 %CU PY 
CALL GAUSEL (A eK sCeDeN) 
OO 8 I=1eN 

8 U(TI=CC1) 
IF (MARK eEQ. 0) GO TO 10 
OIST=V2*DX 
GRAD=(U(N)-U(N1L) )/ZO0IST 
UX=UC(N) 
O00 % T=NeN2 

Ss UCT )=UCN1] )¢GRAD*( I-N1 ) 0X 
UCN241)=UX 

10 NI=N2 
Vil=V2-IREM 
MAR K=0 
OT=DX/(NG*FDERIV) 
TIME 2=TIME/S2.6 
IFCOTeLTeTIME2) GO TO 30 
OT=TIME2 

30 IFC TIMEeLTeTOUTP) GO TO 25 
TOUTP=TOUTP+TOUT 
CALL CUTPUT (UeXeNlo Vi eGAMD peTIME.CV) 
GO TQ 25 
END 


SUBROUTINE GAUSEL (A eG 0CeDsN) 


c----- RE RK EKA MAKE HEE EAE REE EK REKEEESKEKEKE RSE KERR EE EEKRRKEREKE KEKE 
C----- THIS SUBRCUTINE SOLVES *N® LINEAR SIMULTANEOUS TRIDIAGGNAL 
C-----EQUATIONS»s ANO STCRES THE RESULT IN THE ARRAY C(I) cocccee 
Ca HK KEKE KEE AKKEE KEK EKK SERRE EKSEE ESET EK EKEKRKE KE RKKKKKRERKKEEE EKER EKEKERE 


REAL A(1)2841)sCC1)I.001) 
IFCNeGT21)GO TO 8 
CINI=D(NIZA(N) 
GO TO 6 
C----- TO CHANGE TO UPPER TRIANGULAR FORMe 

8 O00 S I=2eN 
ACI) =ACT)—BC1—-1) #CCI)ZACI-1) 
OC1I=OC1T)-DOC 1-1 €CK TI ZACI-1) 

S C(I)=06 


C-=-=—— BACK SUBSTITUTIONee 
C(N41)=O00 
I=N 


4 CIV =(OCT)-8C 1) *CCL4+1)7ZATT) 
IF(LeLEe1)GO TO 6 

7 t=1I-1 
GO TO 4 

6 RETURN 
END 
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SUBROUTINE OUTPUT (UsX »NloV1eGAMDeo TIME eo CV) 
REAL U(1).-UNUT(SO0) 

COMMON DOX.PO 

ALI=Ni-1 

SUM=0.6. 

OO 1 J=2eNL1 

SUM=SUM#22%U(1) 

SUM=SUM4#U(N1L) 

SUM=SUM#*0 .5*DX 

SUM=SUM#0 .S*#VILROXK(UCAL )DFUCNI 41 ) ) 
SE TEND=P0*¥X+GAMD F0eS ¥X¥X 
SET=SUM/SETENC 

SET=1--—StT 

TIMDAY=TIME 

RICAY=SQRT(TIMDAY } 

WRITE (628) XeTINDAYsRTICAY 


P=PO0#GAMD eX 
RETURN 
END 


2/2 


8B FORMAT(#0O*. ® THE DEPTH OF THAW IS %eFil2e4e*°AND T=* oF 1004. ° YEAR 
1S AND ROOT YERRS=" oF1064) 
TIMFACH=CVETIMES(X®X) 
WRITE (6Ge2dITIMFACeSEToNieVI 
2 FORMAT( 40% e* TIME FACTCR CV eT/X¥¥2 =*,F1I4e8,* DEGREE OF SETYTLEMEN 
IT =",F10-5/°*N1I=",15e° AND THE FRACTION [S®* oF1064) 
N=N1 41 
P=P0+GAMD€X 
DO 4 T=1.5N 
4 UOUT(T=UCII7P 
WRITE(6,3) (UGQUT( I) +f=2 oN) 
3 FORMAT(® *©,10F10e¢3) 
RETURN 
END 
SUBROUTINE THAW(XeDERIVeTIME se CBe CN) 
Cc----- THIS SUBROUTINE CALCULATES THE DEPTH OF THAW FROM THE FORMULA 
C----- X(T) = BeT*¥¥*¥N 
C----- THE RCUTINE ALSO CALCULATES THE DERIVATIVE OxX/DTececee 
Cm HH LER EEE KE REK ETEK EEE EEE EKER EK REE KE REE EKER KEKE EEE ERK EEE EKKE EEK 
X=CB*¥(TIME4*CN) 
DER IV=CN¥C3e(TIME*##(CA—1.)) 
RETURN 
END 
SUBROUTINE POSIT( XeoDOXsNeVeNP) 
C= ——— YHIS FOUTINE GIVES THE FOSITICN OF THE FREEZE-THAW INTERFACE 
C—-----F RCM THE RESULTS FROM SUBRCUTINE*THAW*® AND OUTPUTS A NODE AND 
C----- THE FRACTION GF A NODE SPACING WHERE THE INTERFACE ISecccee 
NP=X/DX 
V=X/0X—-NP 
N=AP4+1 
C----- *N*® ITS THE NUMBER OF FINITE OIFFERENCE EQUATIONS FOR THE NEXT 
C----- TIME STEPec cece ececece 
RETURN 
END 
REAL FUNCTION P(PO0+eGAND.X) 
C----— THIS FUNCTION CALCULATES TRE EFFECTIVE OVERBURDEN STRESS 
C----- AT THE FREEZE-THAW INTERFACE ccce 
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The equations for the lower layer are identical except to replace 
By by Bos where 
WD 


B, = 
2 © Tax)? 


These equations apply to the nodes 

i = (m+ 1), ---- (n - 2) 
and the equations for the nodes (n - 1) and n are identical with those 
for the same nodes described in Appendix A.3. 

The only remaining node is the interfacial node 'm', between 
the two layers in question. At this node, the continuity of water flow 
and the equality of excess pore pressures must be satisfied 
Simultaneously with the governing consolidation equation for each layer. 

The node m forms the last node in layer 1 and the first node 
in layer 2. As the excess pore pressure in each layer must be the same 
in both layers at this point, the sets of finite difference equations 
for each layer may be linked by considering these two nodes as the 
same node, and so obtaining only one finite difference equation for the 
node. 

The equations of consolidation for each layer bordering on the 


interface may be written down according to the Crank-Nicholson method 


as follows. For layer 1 


and for layer 2 


uate Umaine Bo (uaa 54172 pe yon Ug j+1/2) (A.25) 


where pandq are two fictional nodes introduced as shown in Fig. A.1. 
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pore 
pressure 


u(x,t) 
- (m- 2) 
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ly LAYER 
(m- |) | 
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pk 
LAYER 
2 
K,1 Cy 
depth 
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Fig. A.1 Finite Difference Grid at the Interface Between 
Two Soil Layers 
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We must now satisfy the continuity of pore fluids at the inter- 
face. In order to be consistent with the Crank-Nicholson difference 
method, central difference operators must be used to approximate the 
first order derivatives in the continuity equation, and this requires 


the presence of the nodes p and q. As central difference operators are 


applied, the discretization error of the Crank-Nicholson method (ax)? 
is preserved. At the time level j+1/2, we have therefore 
1 (in layer 1) = K, 4 (in layer 2) 
m,jtl/2 ~|m, j+1/2 (A.26) 
In discrete form, eq. (A.26) becomes 
most 2 eased ie ayy cred ely2. agit 1/2 
] 2 AX, 2 2 AX, (A.27) 


Equations (A.24), (A.25) and (A.27) now contain three unknown 


values of the excess pore pressure Uy ge1/2? Yq, 4+1/2 and Un jt /2° 


Eliminating Uy j+1/2 and Ug, jt1/2 from these equations, and setting 


ee Mad ) 


Ui, j41/2 igh fated 
for consistency with the Crank-Nicholson scheme, the following implicit 


finite difference equation is obtained for node m: 


~ Uap, ge1 * CTFN/ByAT/BotT A) Un jet 7 Ey Une gH 
= uy? (1/BytE /Bo-I-1) ug * 1 UT (A. 28) 
Karax 
_ Ky Axo 
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Hence the finite difference equations for the two layers are 
linked by another linear tridiagonal equation. The complete set of 
linear equations for the layered profile may now be written in matrix 


form, and solved by the method described in Appendix A.3. 


A.6 Program Listing for Thaw-Consolidation in a Two-Layer Soil Profile 


Purpose:- To determine the excess pore pressure distribution and degree 
of settlement in a thawing soil with a constant applied surface temperature. 
The soil profile is assumed to exhibit two distinct sets of thermal and 


geotechnical properties. 


Description:- The rate of thaw and the resulting thaw-consolidation 
solution are well established for the thawing of the overlying layer, 
and this program determines the solution subsequent to the thaw plane 
entering the underlying layer. The rate of thaw through the lower layer 
is calculated using the simple Stefan-type relationship derived for a 
two layer system in Section 2.9. The program is set up to calculate 

the pore pressure distribution in the two layer system using a Crank- 
Nicholson finite difference procedure similar to that used in Appendix 
A.3. The major difference is that this routine must incorporate extra 
finite difference equations to satisfy the excess pore pressure 


conditions at the interface between the two soil strata. 


Usage:- The program requires data on three cards as follows. 


(i) In FORTRAN format (7F10.0) the necessary parameters in 


order are 
k ee ed | ‘ 
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where kK, and ko are the thermal conductivities of the upper and 
lower layers respectively, 
H is the height of the upper layer, 
T. is the constant applied surface temperature, 
L is the latent heat of the soil, 
te is the termination time for the program, 
Cc, is the consolidation coefficient for the underlying layer, 
Ax is the finite difference interval in the x-direction, 
a and y' are the surface loading and submerged unit weight 
of the soil. 
(ii) The second data card requires the following data in 
order in format (5F10.0), 
H 


Ky, K 


1? 2° Cy? Qs 
where Ky and Ky are the permeabilities of the upper and lower layers, 


( is the consolidation coefficient for the upper layer, 


vl 
and a is the alpha value for the upper layer. 
(iii) The third data card requires the parameter NG in FORTRAN 


format (13), and the description of this parameter is given in Appendix A.4. 


Units:- A suggested set of units is lbs, feet and years. 


: : 2 
Output:- The program prints out the thaw depth, time, time !/ » degree 
of settlement and the normalised excess pore pressure profile, and this 


information is output two hundred times between t = 0; and t = te. 


Example of data input:- For an example of 10 feet of fine-grained soil 


overlying a coarser material, a typical set of data is 
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k 
2 i! 

ery scam] OSEC. 
Cakat 

2 Bom 100 ft/year|/¢ 
Ce = 16 years 
Cy = 225 ft“vyear 
Ax frn5 ft. 
St = 0 
y' = «60 Ibs/ft® 
K, = 1 ft/year 
K, = 4 ft /year 

O 2 

Cy = 25 ft /year 
0 = 10 ft/year!/@ 
H =07 10 ft. 
NG = 50 (See overleaf for program listing. ) 


A.7 Finite Difference Scheme for Thaw-Consolidation with Sand-Drains 
Writing eq. (6.5) and (6.6) in finite difference form according 
to the Crank-Nicholson scheme we obtain respectively 
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c----~ PROGRAM TO CALCULATE THE PGRE PRESSURES IN A THAWING SGILeo 

c----- THE PROGRAM USES THE CRANK-NICHOLSUN FINITE DIFFERENCE METHOD.s. 
C-----IT 1S A FIXED GRID» MOVING BOUNDARY SITUATION» AND THE SCIL IS 
C-----ASSUMED HOMOGENEGUSe.s CNE DIMENSIONAL FLOW IS PERMITTED. AND 
C----- THE PROGRAM CAN INCCRFERATE ANY SPECIFIED RATE OF THAWeeee 

Cc ——— oe 

C-----THE PROGRAM IS NOW ADAPTED TO SOLVE A TWO-LAYER PROBLEM USING A RATE 


C-----OF THAW CALCULATED FRCM THE FORMULA GIVEN BY NIXON €& MCeROBERTS( 1973) 
CHR ROKER EK EE EE CREE EERE EEE KORE EK EE EEEK EER EERE RE KE HEKEKAEE ERE KEKKKES 
REAL U(S00) .A(S00)48( 500) oC (S00) eD( 500) 
CCMMON OX.P0.0H 
C------ ENTER THE PARAMETERS NECESSARY FOR SOLUTION 
99 READ (50101)CB eCNo TF eC VeDX eo FO »GAMD 
101 FORMAT(7F10.0) 
WRITE (60102) CBeCNeTF eCVeDXeP0eGAMD 
102 FORMAT(*#1%,¢ DATA INPUT: —*78 THE RATE OF THAW IS GIVEN BY 
13*/20Xe *P =*,FBe2e° Q= *oFBe3/* THE SOLUTION WILL BE TERMINAT 
2EDO AT TIME=*,FIO0.3/* THE COEFFICIENT OF CONSOLIDATIUN I1S*.F10.5S/ 
3*THE GRID SPACING IN THE X-DIRECTION IS*,F7e3/*°THE APPLIED LOAD AT 
4 THE SURFACE IS*%.eF1064/*THE SUBMERGED UNIT WEIGHT OF THE SOIL IS*, 
SF 8-3) 
REAC(Se35)PK1 oPK2 eC V1 eALPHAH 
35 FORMAT(5F10-0) 
WRITE (6236) PK1 »«PK2eCV1l» ALPHAsH 
36 FORMAT(*O*.® THE PERMEABILITY OF THE UPPER LAYER 15% ,F100e6/% THE PER 
IMEABILITY OF THE LOWER LAYER IS%eF10-.6/*° THE COEFFICIENT GF CONSOLI 
2D0ATION FOR THE OVERLYING LAYER [5% .Fl0¢54% FT*e*2/VYEAR*/*®* THE ALPHA 
3VALUE FOR THE UPPER LAYER IS*% oF 100e50% FISYR¥¥0.5°9/*THE HEIGHT OF 
4THE UPPER LAYER I1S% oF1l00«4% ) 
REAC(S.31)NG 
31 FORMAT(1I3) 
W=PK2/PKI1 
TOUT=TF/200.~ 
TOUTP=TOUT 
M=H/DX 
XO=121*DX 
TO=((X04+CB) *(X0+C8B)—-CB*CB) /CN 
TIME=TO 
OT=TO/20 6 
x=X0O 
DER IV=X/TIME 
CALL POSIT(XeDXeNoVeNP) 


Ni=NP4M 

Vi=V 

N=N141 
Ca KKK EEE ERE KE ER CEH EEKE REE ER EE EEE 
c----- CALCULATE THE INITIAL VALUESeecees 


R=ALPHAZ (22. *SCRT(CV1)) 

RP=SOQRT (32141593) 

AA=PO/ (ERF(R) FEXP(—-R*¥R) SCRP#R)) 
BB=GAMD/(1e 41 e/(2¢4R*R) ) 

00 37 I=1»™M 

Z1i=I 

Z=Z71/N1 


37 UCT J=AA*ERF(R*Z)+BB¥I #0X 
GRAD=PK1*®(U(M)—-UCM—1))/7PK2 


U(M+1)=GRAD4U(M) 
U(M42)=U(Mt1)+4GRAD*V1I 
CALL CUTPUT (UsX oN1 oV1 eGAMD ol IME ) 
25 TIME=TIMEFOT 

IF CTIME*GT oTF) GO TC G9 
BTI=CVL*¥OTS( 0X 0X) 
BT2=Cv*eOT/s(OX4*DX) 
BT=BT1 

G==——— STORE X AND DX/DT AT THE JTH TIME STEP ececece 
XJ=X 
DERJ=DERIV 
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C-~--—CALCULATE X AND DX/DT AT THE (Ut1)TH TIME STEP ececcee 


_—— 


10 


30 


CALL THAW(X»DERIVeTIMEsCBe CN) 

xJ1=X 

DERJI=DERIV 

CALL POSIT(Xe0OXeN2eV2.NP) 

N2=NP4M 

IREM=0 

IF(N2-EQeN1 ) Go TO 1 

IREM=N2—Ni1 

V2=V2tIREM 

MARK=1 

OCL)=2e¢ ¥(Le—-BT)*UCL IFAT #UC2) 

N=N14¢1 

NL 2=N-2 

OO 2 1=2-sNL2 

BT=BT2 

IF CT eLT eM) BT=BT1 

OCT =OT¥ (UC Tel) 4+UCl-1)) #20¥*(12.-BT)*U(T) 
D(C N21 )=BTEVIFUCN—2)7 (1 64¢V1) #¢1 BT )*UCN-1) + BTKUCNI ZO 1 e#V1) 
OO 3 I=1-eNL2 

BT=BT2 

LF (CT otLT eM) 8T=BT1 

ACI )J=2-*( 1-+BT) 

B6(1)=-B8T 

C(t )=-B8T 

Cl1)=0. 

BT=BT2 

A(N-1)=1.+¢8T 

C(N-1)=-B8T¥*¥V2/(1letV2) 
BCN—1)=-BT/(1e+V2) 

CALCULATE THE COEFFICIENTS FOR THE FINITE OIFFERENCE EQUATION 
AT THE BOUNCARY(NTH) NODE coccee 
BJ=CV*¥DT/(VI*V1I €OX*DX ) 
BII=CV*OT/S(V2A¥*V24DX#DX) 

P J=P (FO »GAMDs XJ) 

PJI=P( PO .«GAMD oXJ1 ) 

CJ=V1 *DX*DERI/SCV 

CJI=HV2*DX*DERVJISCV 
AQCNJ=1-45841*C12-¢4¢CJ1) 

C(NI=-BJ1 

BIN) =064 
DON) =(01--BIJ—-BI*F CI) *FUC(NI*BIFUCN—-1) *B4U1*CUL*PIL #BseCuUePy 
ACM) =1etle/BTIL+W/BT2¢+w 

B(M)=-W 

C(M) =—-le 

O(M)=UCM—1) ¢€C1./8T1+Ww/BT2—1e—-W) *€UCM) + WUC MEL) 
CALL GAUSEL (AeBeCoCeN) 

O00 8 {=1.N 

UCT)=C(T) 

IF (MARK EQ 0) GO TO 10 

OIsT=V2*OX 

GRAD=(UC(NI-UCNL)I/SOIST 

UX=UCN) 

00 5S I=N+N2 

UCL Y=UCNI)4GRAD*CI-NI ) *OX 
UC(N2#1) =UX 

N1I=N2 

V1l=V2-IREM 

MARK=0 

OT=DX/(NG*DERIV) 

TIME2=TIME/26¢ 


IF (CT et TeTIME2) GO TO 30 
OT=TINtE!) 
IFC TIMFeLTeTOUTP) GO TO 25 


TOUTP=TOUTP+TOUT 

CALL GUTPUT (Uex eNl eV1l »sGAMDoTIMEsCV) 
GO) TO 2's 

END 
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SUBROUTINE GAUSEL(A eB eC eoDeN) 


CH EKA ERE KE RE RE EKER AKER SEEKS OEK EK GREK HEREC TE REE EKER ERE EKEEEREK KEES KEKE 
C----—- THIS SUBROUTINE SOLVES *“N*® LINEAR SIMULTANEOUS TRIOIAGOGNAL 
C-----EQUATIONSs AND STCRES THE FESULT IN THE ARRAY Cll )eccecce 

Cm EEK K EEK AKKEK EK EE KER ERKERERAEKER EK EEKESE EE EERE KE EERE ERE KER OKKERHEE ERE 


REAL A(1).801),CCL).DC1) 
IF(NeGT.1)GO TO @ 
C(NI=DCNDZACN) 
GO TO 6 
C----—— TU CHANGE TO UPPER TRIANGULAR FORM. 
8 DOS I=2~eN 
ACTI=ACT)-B8C1-1)*CCT) ZACI-1) 
D(1)=OCT)-O¢ 1-1) ¥C( IT) ZACI-1) 
S§uehs0s 


C----- BACK SUBSTITUTICNoe 
C(N41 )=06 
[=N 


4 CCPI=HCOCII-BCLI*CCL +l I/ACT) 
IF(TeLE«t )GU TO 6 


7 t={-1 
GO TO 4 

6 RETURN 
END 


SUBROUTINE OUTPUT (UsX «Nie V1leGAMDseTIMEeCV) 
REAL U(1).U0UT{(500) 
CCMMON DXeP0eH 
NLI=N1i-1 
Z=X+t+H 
SUM=06 
OO 1 T=1,NL1 
1 SUM=SUM#2-¥*U( 1) 
SUM=SUM4U(N1) 
SUM=SUM*0 -5*DX 
SUM=SUM+40 -S4#V1¥*OX¥( UC N11 )FUCN1 41) ) 

SETENO=P0*¥Z+GAMD*¥0 o-S¥2Z*¥Z 

SET=SUM/SETEND 
SET=1--SET : 
TIMDAY=TIME 
RTOAVY=SOQRTCTIMDAY) 
WRITE (6 08)XeTIMDAYVs RTCAY 
8 FORMAT(*0O* o® THE CEPTH GF THAW IS*% o6Fl2e¢4e°%AND T=* oF 100e4s* DAYS 
1 AND ROOT CAYS=* of 104%) 
TIMFACH=CVETIMES(CX¥X) 
WRITE (6 e¢2) TI MFAC eo SETeN1 eo V1 


(4 FORMAT(#O%e*TIME FACTGR CVeT/SX¥¥2 =*,FI4e8,* ODEGREE CF SETTLEMEN 
LT =" .FI0eS/INI=* oI5e? ANDO THE FRACTION I[S*% sF1064) 
N=N1+41 


P=PO0+4+GAMD *Z 
DO 4 I=leN 
4 VOUT(CT)=UCI)/P 
WRITE (653) (UOUT( I)» T=1e)0eN) 
3 FORMAT(® *,10F1003) 
RETURN 
ENO 
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SUBROUTINE THAW(XeDERIVeTINE + CBe CN) 


C-----THIS SUBROUTINE CALCULATES THE DEPTH OF THAW FROM THE FORMULA 
Cees X(T) = be TREN 

C= += THE RCUTINE ALSO CALCULATES THE DERIVATIVE DxX/DT secs 

C------ SEKREKEREKKEEKEKKK EEK KEKEK EK EEK ERKE BK KEKEKEKEEEKKEKEE EK SKK KKRKEKE KEKE 


X=SORT (CB*¥CBtCN#T IME )-CB 

DERI V=CNZ (26 *SQRT(CBE*C8H+CN#TIME )) 
RETURN 

END 


SUBROUTINE POSIT (XeDXsNeVoNP) 


Ce THIS ROUTINE GIVES THE FOSITION GF THE FREEZE-THAW INTERFACE 
C— = ——— FROM THE RESULTS FROM SUBROUTINE*THAW*® AND QUTPUTS A NUDE AND 
C—— THE FRACTION OF A NOOE SPACING WHERE THE INTERFACE [Scecceece 
NP=X/DX 
V=X/DX-NP 
N=NP+1 
GA *N® ITS THE NUMBER OF FINITE OITFFERENCE EQUATIONS FOR THE NEXT 
(CSS TIME STEPecce ceocscce 
RETURN 
END 


REAL FUNCTIGN P(P0+GAMD».X) 
THIS FUNCTION CALCULATES TRE EFFECTIVE OVERBURDEN STRESS 


C----- AT THE FREEZE-THAW INTERFACE cece 
CCMMON CXeGGeH 
P=PO0¢GAMD¥#¥(X+t+hF) 
RETURN 
END 
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The finite difference equations along the thaw line (i = n+1) 
are somewhat complicated, and are derived by combining the finite 
difference expressions for equations (6.5a) and (6.7). The final 


expression is 


2 
ee rls iat teane’ u 
cyt 4 Cy Neh; +s! 
28, 
art Een aeg retary eres Mates 
ar iy EA ital dik led g a Sal Vea 
r 194 i1,j-l act 4 CY, Cy ( Ar3) 16) 


= eal ' 
where P EO OSsauil X 


The complete set of finite difference equations has been 
programmed for solution and a listing of the program follows. The 
magnitude of the time step At is varied automatically to optimise the 
program processing time. The format statement for data input is self- 
explanatory, and the input parameters necessary are,in order; 


1 ' 
)> Y's te 


a, bs Ms, ns Qs Cy? (Re + Go 


where te is the time at which the program will terminate execution. 


A suggested set of units are Ibs., feet and years; or tonnes, metres 


and years. 


A.8 Program Listing for Thaw-Consolidation With Sand Drains 


Purpose:- To determine the excess pore pressures and degree of 
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consolidation settlement in a thawing soil, when the rate of consolidation 


is accelerated by the presence of vertical sand-drains. 


Description:- The program solves a two-dimensional consolidation problem 
in co-ordinates r (radius) and x (depth). The Terzaghi-Rendulic 
equation is solved by an ADI procedure, subject to the boundary condition 
on a thaw interface which moves according to the law: 

Xs a/v t 

The consolidation equation is transformed in such a way as to 
place the Poe arnata eyeten in motion, and the thaw plane is made 
Stationary in the transformed co-ordinate system z = x/X(t). Thus, there 
are a fixed number of finite difference nodes in the z direction spaced 
equally at Az. The degree of consolidation is calculated by vertical 


integration of the excess pore pressures. 


Usage:- The program requires nine parameters in FORTRAN format (2F8.0, 
218, 5F8.0), in the following order 
as Disits N's 0s Cys ie y' and te 
where a is the radius of the sand drain, 
b is the effective radius of the sand drain 
( = 0.525 times the spacing), 
m is the number of finite difference subdivisions 
in the r-direction, 
n is the number of discrete intervals in the z direction, 
a is the constant of proportionality between thaw depth 
and time !/@, 


c is the consolidation coefficient, 
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is the surface applied loading, 
is the submerged unit weight of the soil, 


is the termination time for the program. 


suggested set of units is lbs, feet and years. 


data input: - 
me0.5, ts 
Dac Diehlae 


10 


10 


10 ft/year /¢ 


50. ft“/year 


500 1b/ft® 
= 50 1b/ft° 


16 years 
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REAL A(S0)2E(50)5C{50) 40(50) 0U( 50250) oe V( 50450) 
COMMON MieN1 


1 REAC(S22) ADe Blo MeNo ALF 0 CVeFGe GAMDe TE 
2 FORMAT (2F 8.0021 8,5F 8e 0) 
WRITE (603) ADs 8D eMoeNeoALF oCV eo PO eGAMDe TE 
3 FORMAT( 91 %e *A=* 4F10 03 olLOXs *B="eF I Oo del 0OXe *M=* eo l3olOXe *N=*el3/ 
1° ALPHA=* »Fl1003010Xo *CV=" F106. 34 PPO=* pFIL0e3e10X0® GAMMA SUB o=* pFli0-3 
2/7*°THE FINAL TIME 15%.F102-3) 
OUT=0.-25 
XOUT=OUT 
NG=50 
X=BO/206 
OX=X/N 
TIME=X*€®X/ (ALF ¥ALF ) 
OZ=le/hN 
DR={BCO—-AC)/M 
N1I=N¢1 
MI=M+el 
M2=Me2 
C—-—-—-—-CALCULATE THE INITIAL PCRE PRESSURES cocce 
RR=ALF/(2e*SQRT(CV)) 
SP=SQAhT (32141593) 
AA=P0/(ERF (RRIF EXPL —-RR*RRI/(SP¥RR) ) 
BB=GAMD/(1.+16/(2e¥#RRERR) ) 
O00 @ [=e2.Nl 
KXX=( T-le ) #X/N 
UX=AA#ERF(XX/(2¢*SQRTI(CVETIME))) + BB¥XX 
DO 4 J=leM2 
4& ULT oJ J=UX*( J—-1Le)/SM 
CALL OUTPUT (UeOXeXeTIMEePO0eGAMD ) 
DTI=04eS5/4CCVR(1 of (OX€OX Ft e/(ORFOR))) 
C——---FORM THE EOUNCARY CONDITICNSeecece 
DO 6 I[=1>o6Nl 
U(1-1)=06 
6 V(1eld=065" 
OO 7 J=1 M2 
UC1eJ)=060 
7 V(1leJI=HO60 
5 BY Z=CV*DT/(0Z*0Z) 
BIR=CV*DT/(DR*DR) 
C-----FORM THE LeHe BOUNDARY CONOIT IONe cee 
OO 8 I=1eN1 
8 UCT »M2)=UCLCT mM) 


DO 21 J=1eM1 
21 U(N420J-)= ALF *¥ALF DZ *(PFOtGAMD¥X—U(N1 oJ) I/CVEUINs J) 
DO 10 [=2.N1 
Z=( 1-1.) *0Z 
00 15 J=2eM1 
R={(J—-1e) *ORFAC 
AC J—-1)=1e+2e*BTR 
B(J-1)=-(1.40F/(22*R) 2 *¥BTR 
C(J-1 ) =-(1--DR/(2 o#R) 1 *BTR 
15 D(J—-1)=UC Te J)+ (BIZ/CALF SALFHTIME) )*€ (UCT +1 6 J)—-2e UC ToS) FUCI—-1 0 J) 
14¢0¢25¥4ALF XALFHZHD0Z*(UCT416¢5)-UCT-15I))/7CV) 
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C(M)=—-2.*BTR 
B(M)=0. 


O€C1)I=C§CLI-COL)#*€UC Ll) 
C(1)=0.6 

CALL GAUSEL(A 0B eC oDeM) 
DQ ih J=2eM1 

V(T2»J=C( J—-1) 

CONT I NUE 


—-INCREMENT THE TIME BY ThO FALF TIME STEPSeccece 


TIME=TIME+2.*#0T 
X=ALF*SQRT(TIME) 
CON=8T Z/ (ALF ¥ALFSTIME ) 


C—----ALTERNATE TKE INTEGRATION TO THE Z—DIRECTIGNeecee 


12 


14 


16 
13 


19 


OO 12 I[=2,.Ni 

VCCI oM2)=V(E.M) 

DO 13 J=2.M1 

R=(J—-1.)#OR¢AC 

OO 14 I=2.N12 

Z2=(I-1-)4C2Z 

ACI-1)=1-¢+2e%CON 

B( 1-1 )=—-CON€( 1.40625 ALF €ALF*eZ*0Z/CV) 

C(I -1)=—CONK( 1 e0-Oe25*ALF FALF¥Z*DOZ/CV) 

OCI-1LI=HVC Toe #tBTRECV(LTeJt] —-20*V(C ToS tV( 1s J-2)4+ DR¥e(VOII.Jt1)d 
1-V(CToJ-1)I7(2e#R)) 

ACN) =Le¢2e*¥CONtBTZ¥OZ* (CL etOe. 25 HAL FRALF¥Z*DZ/ CV) /S(CV*ETIME ) 
BC(N)=06 

C(N)=-26*CON 

DIN) =OON) CON ¥#( 1 eo tO e25RALF PALFHEZ*OZ/SCV) *ALF *¥ALF ¥DZ &(POtGAMD EX) SCV 
C(1)=0-6 

CALL GAUSEL (A eB eColoN) 

DO 16 I=2.N1 

UCI esJI=CCI-1) 

CONT {NUE 

OX=X/N 

TF CXeLTeXCUT ) GO TO 19 

xOUT=XOUT+OUT 

CALL OUTPUT(CUsD0XeXsTIMEsPO eGAMD ) 
DI=5eO/(CV¥EC1 oS (COX*OX D1 eo (OR*DR) )) 

IFCTIMEoLTeITF) Go TO § 

co TO 1 

ENO 
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SUBROUTINE OUTPUT (UsDXeXeT IME »P0 eGAMD) 


C----- THIS SUBROUTINE OUTPUTS THE NORMALISED EXCESS PURE PRESSURES 
C~---—THROUGHOUT TFE SOIL MASS» AND COMPUTES THE DEGREE OF SETTLEMENT 
Ce PROFILE ACRCSS A OFAMETER EY VERTICAL LINE INTEGRATIONeoe o 


Ca man — FREER EEK EE HERE EE EEE DEER ARERR TE SHAKE EERE KREE EE EE EERE SE RAKE KEK EKEE EK EEE 
REAL U(S00l) »W(50 250) -SET(S50) 
COMMON Mi1eNI 
WRITE (608)XeTIMEsDX 
8 FORMAT(°O* .* THE DEPTH OF ThAW IS* »pF10e30? AND THE TIME [S*eFl0c4. 
Die OX=".F106«5) 
OO 1 I=1.N1 
P=P0+GAMD¥*¥X 
IF (I 2EQe1) =1e 
OO i J=1leMi 
1 WOT sJ=UC Tes) /P 
WRITE (602) 
FORMAT(*O*».* THE NORMALISED EXCESS PORE PRESSURES ARE3—*) 
DO 3 Y=1eN) 
WRITE (604) (WC Te J) eo J=1 o M1 ) 
FORMAT(*® *,12F10-3) 
CO S J=1oMi 
SUM=(WOC1LeJItWIN1L oJ) IS 20 
N=N1-1 
00 6 I1=2.N 
SUM=SUM4ew (I oJ) 
SET( J)=1e—SUM/N 
WRITE (607) 0SET(C I) eJ=1 oM1) 
7 FORMAT(* *.//* THE SETTLEMENT PROFILE IS* 0/12F10-3) 
RETURN 
ENO 


i) 


WwW 


No 


SUBROUTINE GAUSEL(A eB sC+DeN) 


C----- KRKKEKSERAKERKA KSEE EK TECK KKEKERKKSEK KEK EEKE EEK KKE KKK KEKE RSE KEK KEK KE EEE 
C-----THIS SUBRGUTIANE SCLVES ®#N* LINEAR SIMULTANEOUS TRIOITAGONAL 

C----—— EQUATIONS» AND STCRES TRE RESULT IN THE ARRAY ClLdecvceece 

Cc----- EREKEEKKKKKEK SEERA SKE EKKKEKE SHE EKEKKE KET ERKEHKKE KKK TREC SE KKEKEKHKER 


REAL AC1)28(1)2CC1).D4(1) 
IF(NeGTe1)GC TC @ 
C(N)J=DON) SAN) 
GO TO 6 
C-—--—- TO CHANGE TO UPPER TRIANGULAR FORM. 

8 DO 5 =2eN 
ACLI=ACL)I—-BC I-21) *CCOT)ZACI-21) 
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